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FOREWORD 


This  report  was  prepared  by  the  Los  Angeles  Division 
of  Science  Applications,  Incorporated,  Los  Angeles, 
California  for  the  Vehicle  Dynamics  Division  of  the 
Air  Force  Flight  Dynamics  Laboratory,  Wright-Patterson 
Air  Force  Base,  Ohio.  The  research  was  conducted 
under  Project  1370,  "Dynamic  Problems  in  Flight 
Vehicles,"  Task  137004,  "Design  Analysis,"  Contract 
F3 36 15-74 -C-3094.  James  J.  Olsen  and  later 
Lt.  William  L.  Holman  of  AFFDL/FYS  were  the  Air  Force 
task  engineers. 

R.  M.  Traci  was  the  principal  investigator  for 
the  study.  J.  L.  Farr,  Jr.  contributed  to  the  development 
for  the  numerical  solution  method  and  developed  the 
necessary  computer  programs.  Consultant  E.  D.  Albano 
contributed  to  and  guided  the  research. 

The  authors  submitted  this  report  in  July  1975 
for  publication  as  an  AFFDL  technical  report  to  cover 
research  performed  from  June  1974  through  June  1975. 

Other  reports  prepared  and  submitted  under  the 
aforementioned  contract  are:  AFFDL-TR-74-37,  "Small 

Disturbance  Transonic  Flows  about  Oscillating  Airfoils," 
AFFDL-TR-74-135 , "Computer  Programs  for  Calculating 
Small  Disturbance  Transonic  Flows  about  Oscillating 
Airfoils,"  and  AFFDL-TR-75-103 , "Computer  Programs  for 
Calculating  Small  Disturbance  Transonic  Flows  about 
Oscillating  Planar  Wings." 
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1.0  INTRODUCTION 


A basic  reason  for  studying  unsteady  flows  is  the  pre- 
diction of  the  effect  of  unsteady  aerodynamic  forces  on  a 
flight  vehicle.  The  transonic  speed  regime  is  particularly 
important  in  this  respect  for  the  same  physical  reason  that 
drag  and  other  steady  forces  are  generally  higher  than  in  sub- 
sonic or  supersonic  speed  regimes.  Another  characteristic  of 
transonic  flows  is  the  potentially  large  phase  lag  between 
boundary  motion  and  induced  unsteady  pressures.  These  effects 
tend  to  increase  the  likelihood  of  aeroelastic  instabilities, 
so  that  transonic  speeds  are  most  critical  for  aircraft  flut- 
ter. There  is  a clear  and  present  need  for  accurate  and  effi- 
cient predictive  techniques  for  unsteady  transonic  airloads. 

This  need  motivated  the  investigation  of  the  perturbation  me- 
thod for  transonic  flows  about  oscillating  airfoils  presented 
in  this  and  a previous  report. ^ The  present  report  contains 
a potpourri  of  analyses,  numerical  tests  and  calculated  re- 
sults which  generalize,  verify  and  exercise  the  method  intro- 
duced in  Reference  1.  By  way  of  introduction,  the  following 
description  of  the  nature  of  unsteady  transonic  flows  and  the 
past  and  present  attempts  at  describing  the  inherently  nonlin- 
ear processes  is  presented. 

The  physical  mechanism  underlying  the  build-up  of  pres- 
sure amplitudes  and  phase  lags  has  been  described  by  Landahl2 
by  considering  the  propagation  of  a low-frequency  acoustic  dis- 
turbance. The  upstream- facing  portion  of  the  disturbance  (the 
"receding  wave")  travels  at  a rate  which  is  essentially  the 
difference  between  freestream  speed  and  the  sound  speed.  In  a 
transonic  flow,  this  speed  is  so  slow  that  acoustic  disturbances 
may  interact  and  build  up  to  a finite  amplitude.  In  order  to 
describe  the  phenomena,  the  steady  or  low  frequency  unsteady 
equations  must  necessarily  be  nonlinear,  thereby  providing  the 
essential  difficulty  in  the  analysis  of  transonic  flow.  It 
has  been  pointed  out  by  Rott^,  Landahl^  and  others,  however, 
that  the  unsteady  equations  may  be  linearized  and  uncoupled 
from  the  mean  steady  flow  for  high  frequency  disturbances  near 
Mach  1.  These  so-called  "sonic  theories"  have  been  developed 
for  both  two  and  three  dimensional  wings  by  a number  of  inves- 
tigators-*' ’ • 5 . Many  of  these  classical  studies  of  the  high 
frequency  equation  are  summarized  in  Landahl* s monograph^.  It 
is  often  the  case,  however,  that  the  nonlinear  low  frequency 


regime  is  of  greater  practical  importance  for  aircraft  flutter. 
Fortunately,  a linear  system  can  be  recovered  for  low  frequency 
transonic  flow,  as  demonstrated  by  Landahl^,  by  expanding  the 


steady  flow.  The  linear  equation  remains  quite  formidable 
since  various  coefficients  are  functions  of  the  nonlinear- 
steady  flow  so  that  numerical  solution  procedures,  of  which 
the  present  work  is  one,  are  indicated. 

Analytical  or  semi-analytical  approaches  to  the  solu- 
tion  of  the  unsteady  flow,  based  on  this  linearization  about 
the  steady  flow,  have  recently  been  presented.  Stahara  and ^ 
Spreiter6  and  Isogai7  have  applied  the  "local  3.xnear i za ion 
method,  developed  for  steady,  near  sonic  flows,  to  the  nonuni 
form  unsteady  system.  The  approach  is  promising  but  its  ap- 
lication  to  general  supercritical  flows  is  questionable. 

Of  more  practical  interest  for  general  applications  are 
the  numerical  approaches  which  have  become  highly  developed 
for  steady  flows  in  recent  years.  Two  basic  finite  difference 
techniques  have  been  used;  known  generally  as  time  dependent 
methods  and  relaxation  methods.  The  time  dependent  methods 
solve  the  complete  set  of  conservation  equations  or  the  un- 
steady full  potential  equation  by  forward  differencing  in  time. 
This  approach  has  been  applied  to  the  calculation  of  steady 
flows  as  the  steady  limit  of  the  time  dependent  solution.  The 
method  can  handle  mixed  flows  and  shocks  within  the  accuracy 
of  the  finite  difference  mesh8.  Time  dependent  calculations 
of  flows  about  two  dimensional  airfoils  oscillating  at  super-  9 
critical  Mach  numbers  have  been  reported  by  Magnus  and  Yoshihara 
and  Beam  and  Warming10.  Such  calculations  are  very  important 
as  they  involve  no  inherent  limiting  assumption  but  they  are 
quite  expensive  computer  timewise.  Ballhaus  and  Lomax 
presented  time  accurate  solutions  of  the  low  frequency  unsteady 
small  disturbance  potential  equation.  Recent  applications  o 
the  method1*  indicate  that  it  could  be  a very  efficient  approach 
to  the  oscillating  airfoil  problem. 

Another  promising  technique  for  both  steady  and  unsteady 
flows,  which  must  be  considered  in  a separate  classification, 
is  the  finite  element  method.  The  method  has  been  apE1^ 
the  oscillating  airfoil  problem  by  Bratanow  and  Ecer  f°r  com 
pressible  flow  and  recently  by  Chan  and  Brashears^  for  tran- 
sonic flows. 

The  other  class  of  numerical  techniques,  which  have  received 

wide  application  to  steady  transonic  flows15 » , a£e  the 

ation  methods.  Ehlers17  and  the  present  authors  have  reported 
the  application  of  relaxation  methods  to  unsteady  transonic  flow 

about  harmonically  oscillating  airfoils.  jj*® ntrates 

parently  similar  except  that  the  work  presented  here  concentrates 
in  a consistent  manner,  from  formulation  through  results,  on  the 
low  frequency  regime.  In  the  approach,  the  flow  is  considered 
as  a small  perturbation  on  the  steady  flow  which  results  in  a 
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pair  of  boundary  value  Jg&SZr'S?  .££ 

unsteady  Perturbatlo^.P  ^ll’disturbance  potential  equa- 

?f5S&S'SS-- 

wave  equation  form  is  f°*  Both  steady  and 

ing  upon  the  nature  of  bv^the  mixed  differencing  line 

relaxation^procedur^f irstVintroducedmfor  steady  flows  by 

Murman  and  Cole10. 

Since  the  are 

Kr-SSrna^tSn  2.0 

generalizations  to  the  basic  method  presented  in 

These  generalizations  an  extet^  flous,  wind  tunnel 

lit: 

sented  in  Section  3-0.  along  concerning  the  accuracy, 

tion  4.0  P-=-‘%^  rt  a £ t ™°"hod.  An  attempt  at  veri- 

sssrsr^syt  Effir1 

with^omparisons^o^the^resent  results^to^those^of^an^1' exact" 

numerical  method9  and  64A410  airfoil  oscillat- 

oulated  results  for  a NACA  64AC06  and  64A41U  air 

ing  in  pitch  plunge  and  con  rol  urface  ^tation^  # ^ 

number  range  0.8  to  1.2  and  r airfoil  angle  of  attack 

are  presented  in  «:®'  rf“iefly  studied  and  the  re- 

and  wind  tunnel  Val1  e^f  c . . The  results  for  unsteady 

suits  presented  in  that  section.  The /^ee  degree  of  freedom 

aerodynamic  coefficients  a e d disCussed  in  Section  7.0. 

Finally  &Section  8?S  s^arizes  study  results  and  comments  on 
tte  ci«ent  state  of  development  of  the  method. 
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FORMULATION  AND  SOLUTION  PROCEDURE 


Small  disturbance  theory  is  the  principal  analytical 
tool  for  all  speed  ranges  and  has  become  increasingly  impor- 
tant in  the  transonic  speed  range  in  recent  years.  The  gen- 
eral theory  including  the  unsteady  small  perturbation  approach 
used  in  this  work  is  summarized  in  Section  2.1.  Generaliza- 
tions of  the  basic  method,  described  in  detail  previously  ' 
to  include  higher  frequency  effects  and  an  extension  of  the 
method  to  three  dimensional  planar  flows  are  included  in  the 
summary.  The  required  numerical  solution  methods  for  the 
steady  and  unsteady  systems  are  described  in  Section  2.2 


2.1  Small  Perturbation  Approach  to  Unsteady  Transonic  Flow 
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Figure  1.  Schematic  of  Airfoil  Geometry  and  Transonic  Flowfield 

The  problem  of  interest  is  the  flow  about  an  airfoil 
(two  dimensional)  or  planar  wing  (three  dimensional)  oscil- 
lating with  various  flexible  or  rigid  body  decrees  of  freedom 
in  the  transonic  speed  range.  The  airfoil  geometry,  flowfield 
schematic  and  coordinate  definition  are  given  in  Figure  1 
above.  Rectangular  coordinates  (x,y,z)  are  fixed  to  the  air- 
foil leading  edge  with  origin  at  the  wing  root  and  U,  M , aro 
are  the  freestream  velocity,  Mach  number  and  sound  speed  re- 
spectively. The  airfoil  has  a thickness  ratio  6,  which  is  the 
airfoil  maximum  thickness  divided  by  its  chord  c,  an  angle  of 
attack  a and  a semi-span  b.  The  assumption  is  made  that  6<<1 


r 


and  a is  of  the  same  order  of  magnitude  as  6,  Also,  the  os- 
cillatory motion  of  the  airfoil  is  assumed  to  be  described 
by  a small  non-dimensional  displacement  e<<5  and  a reduced 
frequency  k=wc/U  based  on  airfoil  chord  where  u>  is  the  fre- 
quency of  oscillation. 

Assuming  inviscid,  isentropic  flow,  the  problem  can 
be  reduced  to  the  solution  of  a single  equation  for  a velo- 
city potential  plus  the  tangency  boundary  condition  on  the 
airfoil  surface.  As  is  well  known,  the  derivation  of  a small 
disturbance  theory  for  transonic  flows  requires  a singular 
perturbation  approach.  The  following  scaling  is  thereby 
introduced: 
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and  the  total  potential  is  expanded  about  the  uniform  flow: 

2/3 
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Retainina  all  terms  of  leading  order  in  the  total  potential 
equation  and  boundary  conditions  results  in  the  following 
form  for  the  unsteady  small  disturbance  system. 
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with  boundary  conditions: 
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where  f o is  the  unsteady  airfoil  shape  function  (Equation 
m belo*)  on  the  upper  and  lower  surfaces  respectively,  and 
where  U denotes  a jLp  in  the  enclosed  quantity  between  y = 0 
and  0+  It  is  noted  that  the  airfoil  tangency  boundary  con 
dition' (Equation  4)  and  the  Kutta  condition  (Equation  5)  are 
applied  in  the  small  disturbance  manner  on  y = 0. 

The  system*  of  Equations  3-6  provides  a formu^tion  of 
the  unsteady  airfoil  problem  in  the  nonlinear  domain,  which 
includes  f lowf  ields  with  shocks.  Certain  terms  in  abo_ 

quency  "approximation^  1h?s  approximation  is  in 

ri^?ILf1!l0w?ti^n^e^Urt^hf^ci  on  the  ac- 
curacy  of  the  theory. 

The  low  frequency  approximation  results  in  significant 
numericIJ  advantage  whi  solved  by  a time  ^dent  technique 
as  pointed  out  recently  by  Ballhaus  and  ““ax  - Such  appro 
to  the  solution  of  the  above  system  have  yet  to  be  fully  ex 
ploited  for  the  unsteady  aerodynamics  problem. 

The  approach  presented  herein  for  solving  the 
system  given  above  (Equations  3-6)  is  to  expand  the  pertu  b 
tion  potential  function  in  terms  of  the  unsteady  boundary  dis- 
tance e«l.  From  this  point  on  all  tildas  <M  will  be 
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dropped  with  the  understanding  that  all  variables  are  scaled 
variables.  Harmonic  boundary  disturbances  are  explicitly 
treated : 


f (x,z , t) 


f (x,z) 
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- , > iftt 

+ ef ^ (x, z) e 


(7) 


and  the  perturbation  potential  is  expanded  as  follows: 


iftt 


4»  (x,y,z,t)  = $°  (x,y,z)  + e $'(x,y,z)e  + 


(8) 


Substituting  this  into  the  perturbation  potential  equation  plus 
boundary  conditions  and  combining  terms  results  in  the  follow- 
ing pair  of  boundary  value  problems  for  4>  and  $ respectively. 
(In  the  following  text,  the  superscript  has  been  dropped  from 

4> 1 . ) 


and 


(K-4> 0 ) <J>°  + <f>°  + 4>° 

' ^x' vxx  yyy  yzz 
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System  9 is  recognized  as  the  usual  formulation  for  steady 
transonic  flow  and  system  10  is  the  formulation  for  the  un- 
steady perturbation  thereof.  It  is  noted  that  the  governing 
equation  for  is  linear  but  of  the  same  mixed  elliptic/ 
hyperbolic  type  as  the  steady  solution.  It  is  also  noted  that 
$ is  in  general  complex  thereby  permitting  phase  shifts  be- 
tween field  quantities  and  the  boundary  disturbance.  As  be- 
fore, underlined  terms  in  system  10  are  neglected  for  a con- 
sistent low  frequency  approximation.  Also  for  two  dimensional 
airfoil  sections,  the  z dependence  on  all  quantities  and  the 
4>r?z  terms  in  the  equations  are  neglected. 

The  main  physical  quantities  of  interest  are  the  pressure 
coefficient  and  airfoil  force  coefficients.  The  pressure  co- 
efficient, defined  in  the  usual  manner,  is  given  by: 
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where  the  steady  and  unsteady  scaled  pressure  coefficients  are 
given  to  leading  order  ir  the  small  disturbance  approximation 
by: 

C°  = -2<t>° , C = -2  (4)  + ik(t>)  (12) 

p x p 


The  formulations  of  the  boundary  value  problems  are 
essentially  complete  with  the  exception  of  the  practical  matter 
of  setting  the  boundary  conditions  away  from  the  airfoil,  which 
depends  on  the  particular  problem;  subsonic  or  supersonic  free 
field,  wind  tunnel  wall  etc.  Asymptotic  far  field  solutions 
to  Equations  ]0  have  been  developed  for  two  dimensional  sub- 
sonic or  supersonic  free  air  or  wind  tunnel  flows  and  for  three 
dimensional  subsonic  flow.  These  solutions  are  described  in 
Section  3.0  as  they  apply  to  the  various  generalizations  of  the 
method  and  are  summarized  for  completeness  in  Appendix  A. 


2.2  F inite  Difference  Solution  Method 


The  numerical  solution  procedure  for  the  boundary  value 
problems  for  the  steady  and  unsteady  perturbation  potential, 
is  based  on  the  mixed  differencing,  line  relaxation  procedure 
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developed  by  Murman,  Cole  and  Krupp  ' . They  pointed  out 

the  essential  ingredient  for  the  success  of  relaxation  pro- 
cedures for  the  steady  transonic  potential  equation.  The  key 
to  the  approach  is  to  account  for  the  local  nature  of  the  flow 
{elliptic  in  subsonic  regions,  hyperbolic  is  supersonic  regions) 
in  the  finite  difference  approximation  to  the  governing  equa- 
tions. The  solution  method  used  in  the  present  work  for  the 
steady  perturbation  potential,  4> 0 , is  patterned  after  the  method 
for  general  lifting  airfoils  developed  by  Krupp1  . Extensions 
of  this  method  for  steady  flows,  to  treat  supersonic  freestream 
flows,  wind  tunnel  wall  effects  and  three  dimensional  flows  are 
presented  in  Section  3.0.  The  version  of  this  technique  im- 
plemented for  the  two  dimensional  unsteady  perturbation  1S  °e- 
scribed  here.  Generalizations  of  the  unsteady  solution  method 
are  likewise  presented  in  Section  3.0. 


The  linearity  of  the  governing  equation  for  <J>  would 
suggest  that  many  of  the  powerful  techniques  of  classical 
analysis  could  be  applied.  The  equation  is  seriously  compli- 
cated, however,  by  the  fact  that  various  coefficients  are 
functions  of  the  nonlinear  steady  potential,  for  which  the 
only  general  solution  methods  are  numerical.  In  view  of  tnis 
complication  and  the  success  of  the  mixed  differencing  relax 
ation  procedure  for  the  steady  potential  equation  of  similar 
type,  it  was  decided  to  use  this  same  method  to  develop  un- 
steady solutions. 


Figure  2,  Schematic  of  Numerical  Solution  Domain 
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Consider  the  boundary  value  problem  for  the  unsteady 
transonic  potential,  with  numerical  solution  domain  indicated 
schematically  in  Figure  2.  In  the  numerical  scheme,  a rec- 
tangular mesh  with  general  grid  line  spacing  is  overlayed  on 
the  solution  field  in  physical  x,  y coordinates.  Uneven  grid 
spacing  makes  it  possible  to  concentrate  grid  points  near  the 
airfoil  slit  and  in  regions  where  rapid  changes  in  the  po- 
tential or  its  derivative  (wing  leading  edge,  shocks,  etc.) 
are  expected  and  to  expand  the  grid  in  a regular  manner  out 
to  the  boundaries  of  the  mesh.  The  required  solution  for  <p°, 
which  does  not  depend  on  <f> , is  solved  independently  on  the  same 
or  a more  refined  grid  and  the  resulting  solution  stored  on 
magnetic  tape.  The  converged  solution  so  obtained  is  then  used 
in  the  solution  process  for  the  corresponding  This  approach 
has  the  benefit  that  <f> 0 need  not  be  regenerated  for  each  un- 
steady boundary  disturbance  or  reduced  frequency  of  interest. 

The  local  nature  of  the  equation  at  each  grid  point  is 
determined  by  the  corresponding  value  of  (K-<J>£)j  j = Vi,j  at 
the  same  grid  point.  Then  if  K-<J>x  > 0 (elliptic!  the  x deriv- 
atives of  4>  are  center  differenced  and  if  K-<j>x  < 0 (hyperbolic) 
the  x derivatives  are  backward  differenced.  Near  the  sonic 
line  at  so-called  "parabolic  points,"  the  center  difference 
form  for  the  x derivatives  of  <J>  are  used.  The  finite  differ- 
ence approximation  to  the  unsteady  potential  equation  are  thus 
summarized  as  follows: 
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Where  Vi  vj  a and  the  elliptic  (superscript  E)  and  hyper- 
bolic (superscript  H)  difference  forms  are  centered  or  back- 
ward differenced  in  the  usual  manner. 

Using  the  above  forms,  the  finite  difference  equations 
are  set  up  for  each  column  (x  = constant)  in  the  grid.  Appro- 
priate finite  difference  approximations  are  made  for  the  body 
boundary  condition,  Kutta  condition  and  farfield  boundary  to 
effectively  close  the  system  of  equations.  The  resulting 
sequence  of  linear  algebraic  equations  is  solved  by  Gaussian 
elimination.  After  each  column  is  solved,  it  is  reiaxed  with 
a variable  relaxation  factor,  depending  on  the  local  nature  of 
the  solution.  This  process  is  repeated  for  each  column  in  turn, 
sweeping  the  grid  from  left  to  right  until  the  change  in  <t>  for 
all  grid  points  during  one  grid  sweep  is  less  than  some  arbi- 
trary small  amount.  Iteration  is  also  required  on  the  unknown 
airfoil  circulation  which  is  always  over-relaxed  and  which  is 
updated  along  with  the  farfield  at  regular  intervals.  A grid- 
halving routine  has  been  implemented  with  considerable  improve- 
ment in  the  efficiency  of  the  method.  Complete  details  of  the 
differencing  and  treatment  of  body  boundary  condition,  Kutta 
conditions  and  farfield  boundary  conditions  are  presented  in 
References  1 and  20. 
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3.0 


GENERALIZATIONS  OF  THE  METHOD 


One  of  the  main  goals  of  the  research  reported  here  was 
to  generalize  and  extend  the  method  and  solution  procedure  to 
treat  higher  frequency  effects,  supersonic  freestream  flows, 
wind  tunnel  wall  effects  and  three-dimensional  planar  wing 
flows.  Pertinent  studies  relating  to  these  developments  are 
summarized  in  this  section  and  typical  calculated  results  are 
presented.  In  particular,  the  development  and  implementation 
of  boundary  conditions  required  for  the  treatment  of  supersonic 
freestream  flows,  wind  tunnel  flows  and  three  dimensional  sub- 
sonic flows  are  discussed  in  this  section  and  summarized  in 
Appendix  A. 


3.1  Higher  Order  Frequency  Effects 

The  previous  formulation  of  the  unsteady  flow  problem 
was  restricted  to  low  reduced  frequencies  k ^ 0(S^/J)  and  as  a 
result  did  not  include  the  effect  of  plunge  (or  the  so-called 
"piston  pressure"  component)  since  this  is  a second  order 
(0(k))  effect  at  low  frequencies.  As  summarized  in  Section  2, 
the  formulation  can  be  generalized  to  include  higher  frequency 
effects  and  therefore  explicitly  treat  the  effect  of  plunging 
motion  on  all  unsteady  modes.  Since  all  previous  results  as 
well  as  the  unsteady  results  (for  pitch  and  control  surface  mo- 
tion) presented  in  this  report  used  the  consistent  low  frequen- 
cy formulation,  it  is  important  to  examine  the  effect  of  the 
high  frequency  terms  on  the  results  for  reduced  frequencies  up 

to  k = . 2 . 


The  generalization  to  the  formulation  is  accomplished 
by  retaining  the  ^tt  term  in  the  full  potential  equation  as 
well  as  the  unsteady  terms  in  the  airfoil  boundary  condition 
and  definition  of  pressure  coefficient.  There  is  some  ques- 
tion about  the  validity  of  the  resulting  small  disturbance 
equation  (Equation  3)  for  general  frequencies  but  it  is  ac- 
cepted by  most  authors  as  the  lowest  order  unsteady  small 
disturbance  equation.  Before  proceeding  with  the  linearized 
unsteady  perturbation  analysis  and  the  comparison  of  the  low 
frequency  to  the  generalized  results,  it  is  instructive  to 
examine  the  governing  equation  in  some  detail  to  point  out 
the  essential  difference  between  the  low  frequency  approxima- 
tion and  the  present  more  general  formulation. 
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3.1.1  Characteristic  Analysis  of  Governing  Equation 

The  governing  equation  is  classified  as  a quasi-linear 
equation  of  hyperbolic  type  with  the  non-linearity  concentrated 
in  the  coefficient  (K-tfcx)  of  the  4>xx  term.  This  coefficient 
is  essentially  the  deviation  of  the  local  flow  velocity  from 
sonic-velocity  and  for  the  purposes  of  the  following  analysis 
let  K-<{>x  = V.  For  any  point  (x,y)  in  the  two  dimensional 
physical  plane  the  equation  can  be  linearized  in  the  usual 
manner  and  is  written  as  follows: 


V 4)  + 
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It  is  noted  that  tildes  have  been  dropped  from  all  variables 
with  the  understanding  that  they  are  scaled  variables.  The 
hyperbolic  nature  of  the  equation  as  well  as  the  essence  of 
the  low  frequency  approximation  can  be  made  apparent  by  exam- 
ining the  "Monge  cones"  associated  with  the  equation  at  a 
point.  Introducing  the  following  transformation: 


5 = x 

n =‘-Hv]y 

T = Vt  + X 


the  canonical  form  of  Eq.  (16)  becomes: 


’5?  + ,»nn  ' “ * 


TT 


(18) 


where 


\1 


1 + — V 

1 n 
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The  characteristic  lines  in  the  (x,t)  plane  are  given  by  the 
section  of  the  Monge  cone  in  the  (C,n)  plane  which  from  Eq. 
(18)  are: 


t * v | ii  | t,  = constant 


or  in  x,t  coordinates: 


~ (1  ± n/|u|)x  = 


constant 


Receding  Nave 


Advancing  Wave 


Figure  3.  Characteristics  of  Small  Disturbance  Potential 
Equation 


The  two  families  of  characteristics  at  a subsonic  (V  > 0) 
point  are  shown  in  Figure  3.  They  physically  represent  the 
locus  of  the  wavefronts  f.©f^  the  upstream  facing  or  "receding" 
wave  and  the  downstream  facing  or  advancing  wave.  The  slope 
of  each  curve  is  the  inverse  of  the  wave  speed  in  scaled  coor 
dinates.  It  is  noted  that  in  the  low  frequency  approximation 
P = 1 so  that  the  slope  of  the  downstream  characteristic  is 
0 and  the  advancing  wave  speed  is  infinite.  Thus  this  approx 
imation  effectively  neglects  the  effect  of  the  advancing  wave 


the  physical  assumption  being  that  the  advancing  waves  move 
rapidly  away  from  the  airfoil  and  therefore  have  a negligible 
effect  on  the  flow  as  compared  to  the  receding  wave  which  re- 
mains in  the  neighborhood  of  the  airfoil.  It  is  noted^that 
near  a shock  or  compression  through  sonic  velocity,  V N 0 and 
the  slope  of  the  upstream  characteristic  -1-  10  so  that  receding 
waves  from  downstream  of  the  sonic  point  accumulate  at  the 
sonic  point  and  can  build  up  to  large  amplitude  disturbances. 
This  phenomenon  is  responsible  for  the  inherent  nonlinearity 
of  transonic  flow  as  described  by  T.andahl2. 

As  the  reduced  frequency  increases,  y becomes  appre- 
ciably greater  than  one  so  that  the  slope  of  the  downstream 
characteristic  increases  away  from  zero  and  the  corresponding 
advancing  wave  speed  decreases.  Thus  the  advancing  wave  be- 
comes more  important  as  the  frequency  increases,  and  could  in 
fact  become  equally  important  as  the  receding  wave  for  reduced 
frequencies  of  order  1.  The  generalized  theory  presented  here 
includes  the  receding  wave  and  would  therefore  account  for  such 
higher  frequency  effects. 

One  final  note  before  leaving  the  characteristic  analy- 
sis concerns  the  genesis  of  the  so-called  "sonic  theories". 

It  is  noted  that  for  very  high  frequencies  (y  °°)  and  for  lo- 
cal Mach  number  everywhere  very  close  to  sonic  l-Ma>->-0  or  V-*-0 
the  upstream  and  downstream  characteristics  and  in  fact  the 
Monge  cone  coalesces  to  a single  vertical  line.  In  this  case 
the  upstream  or  downstream  communication  through  the  <*>xx  term 
becomes  negligible  so  that  the  governing  equation  reduces  to 
the  linear  equation: 


♦ 


yy 


(20) 


which  for  harmonic  disturbances  has  a reduced  wave  equation 
form: 


yy 


(21) 


The  equation  is  parabolic  and  can  be  solved  by  classical  meth- 
ods leading  to  the  theories  valid  for  high  frequency  sonic  flow. 
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3.1.2  Unsteady  Perturbation  Analysis 


The  development  of  the  linearized  unsteady  perturba- 
tion for  the  general  unsteady  small  disturbance  system  pro- 
ceeds in  the  same  manner  as  described  in  detail  for  the  low 
frequency  system  in  previous  work!  and  as  summarized  in  Sec- 
tion 2.  The  resulting  system  of  equations  (Eq.  10)  varies 
from  the  low  frequency  theory  by  way  of  the  addition  of  fre- 
quency dependent  terms  in  the  governing  equation,  body  bound- 
ary condition  and  wake  condition.  The  treatment  of  these 
terms  in  the  finite  difference  solution  method  is  straight- 
forward and  they  have  been  added  to  the  computer  programs 
developed  in  the  study  as  an  option  thereby  facilitating  the 
comparison  of  the  low  frequency  formulation  to  the  present 
formulation.  Some  comments  are  contained  in  this  section 
concerning  the  implementation  of  these  various  changes  to  the 
method . 


The  governing  equation  for  the  unsteady  perturbation 
potential  requires  but  the  addition  of  the  kfi<}>  term  which  in- 
cludes the  advancing  wave  in  the  formulation.  The  term  is 
included  in  the  finite  difference  equations  in  the  obvious 
manner.  Of  greater  practical  importance  to  the  results  of 
airfoil  pressures  at  low  frequencies  are  the  addition  of  un- 
steady terms  to  the  body  boundary  and  wake  condition  and  to 
the  definition  for  pressure  coefficients.  This  is  because 
these  terms  directly  affect  the  flow  near  the  airfoil  whereas 
the  modification  to  the  field  equation  indirectly  affects  the 
flow  about  the  airfoil.  For  example,  the  new  body  boundary 
condition  is: 


<P 


y 


+ ikf 


(22) 


It  is  noted  that  the  low  frequency  approximation  neglects  the 
effect  of  boundary  motion  on  the  downwash  including  only  the 
effect  of  the  quasi-steady  change  in  airfoil  slope.  As  noted 
earlier,  the  pure  plunging  mode  (fE  = -1)  results  in  the  null 
solution  ( = 0 ) for  the  low  frequency  formulation  but  is  in- 
cluded in  the  present  formulation.  It  is  expected  that  the 
addition  of  the  plunge  motion  to  the  downwash  at  the  wing  will 
have  its  greatest  effect  on  the  out  of  phase  component  { Xm<J> ) 
of  the  solution  and  should  therefore  affect  the  phase  of  the 
resulting  airfoil  pressures.  Based  on  the  general  form  of  the 
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body  boundary  conditions  and  definition  of  pressure  coeffi- 
cient (below)  it  is  expected  that  this  effect  will  be  of 
0 (k)  . 

The  other  change  in  the  formulation  involves  addition 
of  the  unsteady  term  to  the  equation  for  the  small  disturbance 
pressure  and  the  resulting  effect  that  this  has  on  the  wake 
condition.  The  equation  for  the  unsteady  component  of  the 
scaled  pressure  coefficient  is: 


C = -2(4>  + ik(J>)  (23) 

P X 


Again,  for  low  frequencies,  the  primary  effect  of  the  general 
ized  formulation  will  be  on  the  out  of  phase  (Im  Cp)  component 
of  the  pressure. 

The  Kutta  condition  requires  that  the  jump  in  pressure 
at,  and  downstream  of,  the  airfoil  trailing  edge  be  zero: 


[<f>x  + iktfjg*  = A4>x  + ik  A<)>  = 0 (24) 

on  y = 0,  x Z 1 


where  A4>  is  the  jump  in  potential  across  the  wake.  Eq.  (24) 
can  be  solved  explicitly  to  give: 


A<t>  = A4>te 


ik(x-l) 

e 


(25) 


y = 0,  x £ 1 


where  A4>te  is  the  jumP  potential  at  the  trailing  edge  equal 
to  the  airfoil  circulation.  It  is  recalled  that  in  the  low 
frequency  approximation  At})  = A<j>te  s°  that  this  generalization 
involves  an  oscillatory  change  in  phase  of  the  wake  jump  con- 
dition due  to  the  downstream  facing  wave. 


— — - 
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The  generalized  forms  for  the  body  boundary  condition 
(Eg.  22)  wake  condition  (Eq.  25)  and  pressure  coefficient 
(Eq.  23)  are  incorporated  into  the  finite  difference  equa- 
tions in  an  analogous  fashion  to  the  low  frequency  version 
discussed  in  detail  in  Reference  1 and  they  need  not  be  re- 
peated here.  The  finite  difference  solution  method  also 
proceeds  in  the  same  manner  and  the  comparison  calculations 
presented  in  the  next  section  indicate  that  the  generaliza- 
tions described  here  have  little  if  any  effect  on  convergence 
or  other  numerical  details  of  the  scheme. 

3.1.3  Results  for  the  Generalized  Formulation 

Results  are  presented  in  this  report  using  the  low 
frequency  approximation  for  various  unsteady  modes  of  oscil 
lation  and  for  reduced  frequencies  up  to  k = .2.  It  is  ap- 
propriate to  verify  the  accuracy  of  these  results  in  light 
of  the  higher  order  frequency  effects  discussed  above.  To 
this  end,  calculations  have  been  performed  for  a flat  plate 
oscillating  in  pitch  in  a subsonic  stream  and  for  a NACA  64A410 
airfoil  oscillating  in  pitch  and  plunge  at  a supercritical  Mach 
number.  Results  for  the  low  frequency  formulation  are  compared 
to  the  modified  formulation  and  to  results  of  alternate  formu- 
lations . 

Results  for  the  flat  plate  airfoil  oscillating  in  pitch 
[ (f '=-l-ik(x-xp)  ] at  MOJ= . 7 and  a reduced  frequency  k = .2  are 
presented  in  Figure  4.  Distributions  of  the  Real  and  Imaginary 
parts  of  the  jump  in  pressure  coefficient  across  the  airfoil 
per  unit  angle  of  attack  are  shown.  The  results  calculated  by 
the  present  method  in  its  approximate  low  frequency  and  gener- 
alized form  are  compared  to  the  corresponding  results  calcu- 
lated by  the  doublet-lattice  method21.  The  calculations  are 
for  perturbations  about  a uniform  (<t>°  e 0)  mean  flow  and  the 
generalized  results  should  compare  within  numerical  error  to 
the  doublet-lattice  results.  Indeed,  as  shown  in  the  figure, 
the  results  of  the  high  frequency  formulation  compare  within 
a few  percent  to  the  doublet-lattice  solution.  It  should  e 
pointed  out  that  the  present  method  used  26  finite  difference 
points  on  the  airfoil  and  the  doublet-lattice  method  used  25 
chordwise  boxes  so  that  they  are  being  compared  on  the  same 
basis  in  terms  of  numerical  definition.  Results  using  the 
low  frequency  approximation  are  also  shown  in  Figure  4 and 
they  indicate  the  expected  effect  of  the  high  frequency  terms 
on  the  out  of  phase  component  (Imaginary  part)  of  the  solution. 

The  corresponding  results  for  a NACA  64A410  airfoil 
with  a mean  angle  of  attack  of  2 oscillating  in  pitch  about 
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the  leading  edge  at  k = .2  are  shown  in  Figures  5 and  6.  Fig- 
ure 5 presents  the  unsteady  component  of  lift  as  a function  o'- 
angle  of  attack  (or  time)  and  Figure  6 gives  the  perturbed 
pressure  coefficient  distribution  over  the  airfoil.  In  Fig- 
ure 5 the  results  of  the  low  frequency  and  generalized  fre- 
quency formulation  are  compared  to  the  "exact"  numerical  re” 
suits  of  Magnus  and  Yoshihara9  (see  Section  5 for  more  details) 
Mote  that  for  an  integrated  quantity  such  as  lift  the  higher- 
order  frequency  effects  produce  but  a slight  change  in  magni- 
tude and  phase.  It  is  noted  that  both  formulations  compare 
surprisingly  well  with  the  Magnus  and  Yoshihara  results.  The 
comparison  of  pressure  distributions  in  Figure  6 shows  the 
slight  effect  of  the  high  frequency  terms  on  the  in-phase  com- 
ponent and  a somewhat  larger  effect  on  the  out  of  phase  compo- 
nent; similar  to  the  result  for  the  fully  subsonic  flow.  As 
an  aside,  a comparison  of  Figures  4 and  6 show  the  large  tran 
sonic  flow  effect  for  an  airfoil  with  a practical  thickness, 
camber  and  mean  angle  of  attack.  The  comparison  with  the 
Magnus  and  Yoshihara  results  given  here  and  in  Section  o show 
that  the  present  theory,  in  either  its  low  frequency  approxi- 
mate form  or  its  general  frequency  form,  model  these  non- 
uniform  flow  effects.  These  strong  transonic  flow  effects  are 
likewise  shown  in  Figure  7 which  presents  results  calculated 
using  the  present  method  for  the  same  NACA  64A410  oscillating 
in  uniform  plunge  at  k = .2.  A comparison  of  the  perturbation 
pressures  due  to  pitch  and  plunge  in  Figures  6 and  7,  respec- 
tively show  that  for  low  frequencies,  pressures  and  forces  due 
to  plunge  are,  in  magnitude,  approximately  20  percent  of  the 
pressure  and  forces  due  to  pitch.  This  is  consistent  with 
the  estimate  of  higher  order  frequency  effects  which  indicate 
that  such  effects  are  0(k) . 

In  summary,  the  generalizations  of  the  low  frequency 
theory  to  include  higher  frequency  effects,  necessary  for 
the  treatment  of  plunge  oscillations,  would  seem  to  give  an 
accurate  representation  of  such  effects.  The  comparisons 
presented  here  show  that  these  effects  are  0(k)  in  the  low 
frequency  regime  so  that  the  generalized  formulation  should 
be  used  for  all  unsteady  modes  for  reduced  frequencies  as 
large  as  k = . 2 . Since  the  high  frequency  terms  do  not  seem 
to  degrade  the  efficiency,  convergence  or  stability  of  the 
numerical  solution  method,  the  generalized  formulation  pre- 
sented here  would  seem  to  be  recommended  over  the  previous 
low  frequency  formulation. 
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Supersonic  Freestream  Calculations 


The  theory  and  numerical  solution  method  described  in 
the  report  applies  equally  well  to  subsonic  or  supersonic 
freestream  flow.  The  application  of  the  method  to  supersonic 
flows,  the  subject  of  this  section,  requires  that  certain  de- 
tails of  the  method,  in  particular  the  treatment  of  the  far- 
field  condition,  be  modified.  This  generalization  has  been 
accomplished  for  small  disturbance  steady  flows  by  Murman22'23 
and  his  approach  is  followed  in  the  present  work  as  summarized 
in  Section  3.2.1.  Other  details  concerning  the  treatment  of 
steady  supersonic  flows  which  were  considered  during  the  pre- 
sent research  are  also  discussed  in  that  section.  The  ana- 
logous treatment  of  unsteady  supersonic  flows,  the  prime  ob- 
ject of  this  research,  is  discussed  in  Section  3.2.2. 


3.2.1  Steady  Flows  With  Supersonic  Freestream 


The  primary  generalization  to  the  numerical  solution 
method  for  treating  supersonic  freestream  flows  involves  the 
derivation  of  a farfield  boundary  condition.  The  asymptotic 
solution  available  for  subsonic  freestream  flows  (Appen- 
dix A)  clearly  do  not  apply.  Even  if  such  a solution  were 
convenient  to  derive  for  supersonic  flows,  the  application 
of  a Dirichlet  boundary  condition  is  not  appropriate  and 
could  lead  to  numerical  difficulties  for  a hyperbolic  far- 
field flow.  The  method  used  here  involves  developing  asymp- 
totic approximations  to  the  characteristic  relations  for  the 
hyperbolic  form  of  the  transonic  potential  equation  (Eq.  9)  , 
and  to  apply  these  relations  as  boundary  conditions  on  the 
far  boundary  of  the  grid. 

Letting  u = and  v = <f>y  the  characteristic  directions 
and  relations  for  the  two  dimensional  steady  transonic  poten- 
tial equation  are  respectively: 


i * .(u-K)-1*  (26) 

The  characteristic  relations  can  be  immediately  integrated 
to  give: 


2 3/2 
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where  R"  are  the  Riemann  invariants.  It  is  noted  that  far- 
upstream  the  flow  is  uniform  so  that  u=v=4>H0  as  x-*— °°.  This 
uniform  flow  condition  provides  the  needed  upstream  boundary 
condition  which  is  applied  as  a Dirichlet.  condition  on  4> . 

In  the  transverse  direction  (y -*±°°)  away  from  the  airfoil  the 
characteristic  relations  may  be  approximated  under  the  as- 
sumption that  u<<-K.  Expanding  equation  27  for  small  u gives: 


' tv  + < /-K  u + — — + . 

( 4/-K 
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It  is  noted  that  K<0  for  supersonic  freestream  flow.  For 
y-> ±oo  the  assumption  is  made  that  there  is  no  disturbance  at 
upstream  infinity  (R~0)  so  that  taking  the  first  order  ap- 
proximation to  the  incoming  characteristic  relation  gives: 


4°  = t /^K  cf)0  as  y-±«  (29) 

Thus  as  long  as  the  transverse  boundaries  of  the  finite  dif- 
ference grid  are  supersonic  and  far  away  from  the  airfoil 
Equations  29  provide  an  accurate  representation  of  the  far- 
field  solution. 

This  boundary  condition  is  incorporated  into  the 
finite  difference  procedure  in  the  same  manner  as  the  airfoil 
boundary  condition.  That  is,  on  the  top  and  bottom  of  the 
rectangular  grid  (4>yy)i,j  is  differenced  as  a one-sided  dif- 
ference using  Equation  29  in  the  finite  difference  equation. 
For  example  on  the  top  boundary  (y=ymax) * 


2 

Ayj-1 


for  j=JM 


(30) 


With  4>y  from  equation  29  and  the  backward  difference  form 
for  <t>^Y  (since  the  flow  is  supersonic)  : 
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Using  this  relation  in  the  finite  difference  analog  to  the 
partial  differential  equation  (Eq.  14  above)  effectively 
closes  the  system  of  equations  for  a collumn  of  grid  points 
and  the  solution  proceeds  as  before.  The  downstream  (x-+~)  bound 
ary  of  the  grid  requires  no  special  treatment,  as  long  as  the 
flow  there  is  supersonic,  because  of  the  lack  of  upstream  in- 
fluence in  supersonic  flow. 

In  order  to  illustrate  the  method,  fully  supersonic 
results  are  presented  in  Figure  8 and  9 for  a 6 percent  thick 
circular  arc  airfoil  at  = 1.5  and  a 0.5"  wedge  at  M,„  = 2.0 
respectively.  The  circular  arc  solution  is  compared  to  linear 
supersonic  theory  which  would  be  quite  accurate  for  this  case. 
The  calculation  verifies  the  application  of  the  method  to 
supersonic  freestream  flows  and  demonstrates  that  the  tran- 
sonic theory  reduces  to  linearized  supersonic  theory  for  «>. 

The  present  solution  shows  a slight  degree  of  nonlinearity  in 
this  case  as  well  as  demonstrating  the  numerical  smoothing  of 
the  shock  waves  at  the  leading  and  trailing  edge  of  the  air- 
foil. The  wedge  solution  (Figure  9)  was  performed  to  verity 
the  accuracy  of  the  scheme  for  very  weak  oblique  shocks.  Re- 
sults for  three  grid  designs  with  different  values  of  Ay  are 
shown  which  indicate  that  the  sharpness  of  the  shock  is  strongly 
dependent  on  grid  design.  The  "overshoot"  in  the  shock  Dumps 
is  most  likely  due  to  the  extrapolation  procedure  used  for  de 
fining  values  of  potential  on  the  airfoil  since  it  is  no.  evi- 
dent in  the  shock  jumps  away  from  the  airfoil.  The  overshoo 
is  also  dependent  on  the  strength  of  the  shock  as  shown  bi  th_ 
result  for  a stronger  oblique  shock  in  Figure  10.  Fo^  J-JJ-J;® 
case  with  an  order  of  magnitude  stronger  shock  no  overshoot 
is  evident.  It  is  also  noted  that  the  overshoot  can  be  large  y 
elliminated  even  for  the  weak  shock  case  by  the  proper  choice 
of  Ay  in  the  first  grid  point  above  the  airfoil.  This  ind 
cates  that  a more  intelligent  extrapolation  method  for defin 
ing  the  solution  on  the  body,  such  as,  for  example,  extrapo 
tion  along  characterisitics  (Eq.  26),  would  also  ellimmate 

the  overshoot. 
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To  examine  the  accuracy  of  the  basic  hyperbolic  dif- 
ferencing scheme,  a general  higher  order  finite  difference 
form  has  been  implemented  on  a trial  basis.  The  method  uses 
an  additional  upstream  point  for  the  finite  difference  forms 
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For  v=l  the  scheme  reduces  to  that  presently  used  and  for 
v=0  it  is  a second  order  accurate  hyperbolic  scheme  after 
Richtmeyer  and  Morton23,  Results  for  the  2°  wedge  with 
1.25  are  shown  for  v=0  and  v=l  in  Figure  10  for  comparison. 
Both  methods  give  the  same  shock  jump  and  differ  only  in  the 
sharpness  of  the  shock.  This  comparison  indicates  that  the 
three  point  scheme  used  in  this  work  is  quite  accurate  and 
in  fact  gives  sharper  shock  waves  than  higher  order  accurate 
difference  schemes.  Figure  10  also  shows  the  solution  away 
from  the  surface  (at  y = .2  and  .4)  which  indicates  that  ob- 
lique shock  waves  away  from  the  surface  undergo  considerable 
numerical  smoothing.  Note,  however,  that  the  shock  jump  con- 
dition is  retained.  The  degree  of  shock  smearing  indicated 
is  an  attendent  result  of  any  shock  capturing  numerical  method 
and  can  probably  be  elliminated  only  by  using  shock  fitting; 
a rather  complicaLed  prospect  for  general  transonic  flows. 
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As  described  in  previous  work'*’  and  summarized  in 
Section  2,  the  numerical  solution  method  takes  no  special 
alooin?  of  shock  waves.  As  shown  in  the  above  results  the 
procedure  is  adequate  for  oblique  shocks,  but  as  recent ]y 
pointed  out  by  Murman24,  the  treatment  of  strong  shocks 
(supersonic  to  subsonic  transition)  is  inaccurate.  The  ac 
curacy  of  the  method  can  however  be  improved  by  using  a 
special  "shock  point  operator"  (SPO)  as  suggested  by  Murman 
or  by  using  shock  fitting,  as  implemented  recently  by  Cheng 
and  Hafez25.  since  strong  shocks  are  important  for  transonic 
flows,  especially  low  supersonic  flows,  it  was  decided  to 
implement  a form  of  Murman 's  shock  point  operator. 

A shock  point  in  the  finite  difference  approximation 
to  the  transonic  potential  equation  is  defined  as  a point 
U,j)when  j<0  and  <K-*°)i ,1>°-  such  a point, 

the  finite  dif ference  equation  is  generalized  to. 


0(  K-4>°  )?  ° C + 0(K-«D®)J  + ^yy.  . 

v ^x'lo^xx^j  X 1,3  xxi,j  yyi,D 


where  superscripts  b and  c indicate  backward  and  central 
differencinq.  Murman  suggests  a* 0=1  as  the  form  which  pro 
vides  a "conservative"  finite  difference  scheme . However , 
as  can  easily  be  seen  (and  as  noted  by  Mur^n)  this  is  in- 
consistent with  the  governing  equation  to  O(Ax).  In  the 
present  work  itv*s  decided  to  use  o=l/2,  B-l/2  which  retains 
the  consistency  of  the  finite  difference  approximation  to 
t£s  transonic  potential  equation  but  is  nonconservative  by 
Murman's  definition  to  O(Ax).  This  form  was  chosen  because 
consistency  is  believed  to  be  of  utmost  importance  and  also 
since  the  conservation  property  of  the  finite 
scheme  depends  upon  the  definition  of  the  vari 
quantity  over  the  sides  of  a finite  difference  cell.  This 
can  clearly  only  be  defined  to  first  order  especially  near 
a shock.  In  short  it  is  felt  that  the  shock  pomt 
encing  used  is  consistent  with  the  : formalism  of  shock  cap- 
turing and  could  only  be  improved  by  shock  fitting. 

Results  using  the  shock  point  formalism  are  shown  in 
Fi cures  11.12,  and  13  at  the  end  of  the  section.  Figure  11 
presents  the  pressure  coefficient  along  the  axis  for  a Cir- 
cular arc  airfoil  at  M =1.15.  Results  calculated  with 
and  without  the  SPO  are  presented  and  show  the  imP°r^Jn^  Per- 
fects on  the  shock  jump  and  stand-off  distance.  It  is  i 
estinq  to  note  that  the  pressure  coefficient  on  the  airfoil 
(x/c>0)  is  curiously  unaffected  by  the  treatment  of  the  sh 
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Results  calculated  with  ot— 1,  3—1  (as  per  Murman)  are  indis- 
tinguishable from  the  SPO  result  shown  in  the  figure.  It  is 
noted  that  the  use  of  the  SPO  also  seems  to  have  a beneficial 
effect  on  the  convergence  of  the  iteration  procedure. 

Figure  12,  and  13  are  presented  to  demonstrate  the 
application  of  the  final  developed  method  for  supersonic 
freestream  flows  to  practical  airfoils  of  interest  in  the 
present  work.  The  figures  give  the  airfoil  pressure  dis- 
tribution for  a sharp  leading  edge,  6 percent  thick  circular 
arc  airfoil  at  M^l.25  and  1.15  and  a blunt  leading  edge 
NACA64A006  at  Moo=1.5  and  1.2  respectively.  The  figures  show 
the  expected  non-linear  transonic  results  and  in  the  case  of 
the  NACA64A006  (Figure  13)  shows  the  development  of  and  move- 
ment away  from  the  airfoil  of  the  detached  bow  shock  as  free- 
stream Mach  number  approaches  Mach  1. 


Airfoil  Pressure 


Figure  9.  Pressure  coefficient  for  Wedge  in  Supersonic  Flow 


Figure  10.  Pressure  Coefficients  for  a Wedge  in  Supersonic  Flow 


Figure  12,  Airfoil  Pressure  Coefficients 


3.2.2  Unsteady  Flows  With  a Supersonic  Preestream 

The  extension  of  the  unsteady  perturbation  method  to 
flows  with  supersonic  freestreams  proceeds  in  an  analogoas 
manner  to  the  steady  flow  case  just  discussed.  The  main 
generalization  involves  the  development  and  implementation 
of  a farfield  boundary  condition  for  the  numerical  solution 
scheme.  As  in  the  steady  case,  a characteristic  condition 
is  developed  which  is  valid  far  from  the  airfoil.  To  this 
end,  consider  the  governing  equation  for  the  unsteady  per- 
turbation potential  written  in  the  following  form: 


4> 


xx 


— - <j> 

I K | ™ 


(21a*++i*x) 


(35) 


where  it  is  noted  that  K<0  for  supersonic  freestream  flow. 
Introducing  the  characteristic  coordinates: 


C = x + S\K\  y 


ri  = X - >/ 1 K | y 


(36) 


Equation  (35)  can  be  written  in  the  following  canonical  form: 


3 z(j)  _ 1 

3£3n  " 4~[~kT 


(37) 


Integrating  along  the  £ and  n characteristics  in  turn  results 
in  the  characteristic  relations: 


+ (2ift<|>+<|>£<j>x)  dn} 


(38) 


dy  _ l 

on  - -jm 
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and 


3<J> 

3n 


- /{2in*+«®*  ) 
4 | K | < X X 


+ 


(2ifl<}>+4>  °<f>  ) d£ 
x A n 


(39) 


on 


dy  1 


As  before  upstream  infinity  x-*— ■ 00  is  considered  a region  of 
uniform  flow  so  that  $«$x=0  and  this  is  used  to  set  an  up- 
stream boundary  condition  on  <f>  . On  the  top  and  bottom  of 
the  grid  the  characteristic  condition  valid  along  incoming 
characteristics  given  by  Equations  38  and  39  are  used.  Under 
the  assumption  that  these  boundaries  are  far  from  the  airfoil 
and  using  the  fact  that  the  incoming  characteristic  originates 
in  a region  of  uniform  flow,  the  integrals  in  the  character- 
istic equations  may  be  neglected  so  that: 


<f>  = +Mk!  + — h=  (2i«4>+*°<f>  ) 

Y ' X 2sT\K\  X X 


as  y+±°° 


(40) 


These  boundary  conditions  are  incorporated  into  the 
finite  difference  procedure  by  using  a one  sided  difference 
form  for  <t>yy  at  the  top  and  bottom  of  the  grid  with  the  re- 
quired value  of  <f>y  given  by  Equation  40.  For  example  at  the 
top  boundary: 
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4>  VVi  -i 


yyi,j  Ayj  _x 


- /IkT  ^ - -H=(2in4>  .+4>°b  4>*  ) 

xi,j  2v/j7|'  1'3  Xi,j  Xi,j/ 


(41) 


Ayj_x  (^i , j~^i , j-l)j 


As  before  the  x derivatives  of  <f>°  and  <p  are  backward  differ- 
enced since  the  boundary  points  are  required  to  be  supersonic 
for  this  boundary  condition  to  apply.  The  use  of  Equation  41 
and  an  analogous  equation  for  the  bottom  boundary  effectively 
closes  the  system  of  equations  for  a column  of  grid  points. 

Figures  14  and  15  below  are  presented  to  illustrate 
the  application  of  the  method  to  fully  supersonic  flows  and 
low  supersonic  flow  (mixed  flow)  respectively.  Figure  14 
presents  results  for  change  in  pressure  coefficient  per  unit 
quasi-steady  (k=0)  deflection  in  pitch  (fg=-l)  for  a flat 
plate  (uniform  mean  flow  (j>°E0)  and  for  a circular  arc  air- 
foil both  at  Moo=l ,5.  The  flat  plate  result  should  compare 
identically  with  supersonic  linear  theory  which  is  also  shown 
in  the  figure.  As  shown,  the  present  flat  plate  result  com- 
pares very  well  with  linear  theory  with  the  exception  of  some 
smoothing  of  the  centered  expansion  wave  at  the  leading  edge 
and  the  shock  wave  at  the  trailing  edge.  The  present  theory 
also  shows  a tendency  to  overexpand  the  flow  at  the  leading 
edge  which  is  believed  due  to  the  extrapolation  procedure 
discussed  above  for  the  steady  flow.  The  circular  arc  results 
(Figure  14)  and  results  for  a NACA64A006  airfoil  (Figure  15) 
are  calculated  as  quasi-steady  (k=0)  pitch  perturbations  to 
the  non-uniform  steady  solutions  presented  in  Section  3.2.1. 
The  circular  arc  results  show  an  appreciable  non-uniform  flow 
effect  even  for  a fully  supersonic  flow  of  Mro=1.5.  As  ex- 
pected, this  effect  is  shown  to  be  even  greater  in  Figure  15 
for  the  NACA  airfoil  at  a low  supersonic  Mach  number  (M<»=1.2) 
In  this  case  it  is  recalled  that  the  steady  flow  is  subsonic 
near  the  leading  edge  which  is  responsible  for  the  leading 
edge  singularity  in  the  perturbation  pressure  coefficient 
shown  in  Figure  15. 

The  results  presented  here  demonstrate  the  applica- 
tion of  the  present  method  to  supersonic  flow  and  verify  that 
the  method  reduces  to  linear  supersonic  theory  for  M0O>>l. 

The  method  also  seems  to  give  a reasonable  description  of  the 
transonic  non-uniform  flow  effects  and  this  is  amplified  upon 
in  later  sections  of  the  report. 
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Pressure  Coefficient  Due  to  Pitch  Deflection 


Pressure  Coefficient  Due  to  Pitch  Deflection 


j 


3.3.  Wind  Tunnel  Wall  Effects 

Wind  tunnel  wall  interference  effects  are  extremely 
important  for  steady  transonic  flows  and  because  of  possible 
resonance  effects  are  potentially  even  more  u" 

steady  flows.  The  inclusion  of  such  effects t in  any  theory 
such  as  the  present  one  clearly  enhances  the  t 

theory  for  purposes  of  comparing  with  experimental  data.  Re 
cent  attempts  have  been  made  by  Murman28,  Kacprzynski2?  and 
Collins  and  Krupp28  to  evaluate  wind  tunnel  wall  effects  for 
steady  transonic  flows,  and  an  initial  attempt  at  including 
such  effects  in  the  present  steady  and  unsteady  theories  1 
described  in  this  section. 

plenum  chamber 

. — v,  y,  y , x i x 


i*=IM 


H-  c — H 


■i-1  porous  or  solid  wall  I 

Figure  16,  Schematic  of  Wind  Tunnel  Configuration 

A general  two-dimensional  wind  tunnel  test  configura- 
tion is  sketched  in  Figure  16  above.  It  is  ■ 
airfoil  is  on  the  tunnel  centerline  and  ■ that  1?  The 

porous  wind  tunnel  walls  are  separated  by  a _™e 

treatment  of  such  a configuration  by  the  present cln- 
method  requires  the  inclusion  of  a proper  wall  boundary  con 
dition  in  the  finite  difference  mesh  at  the  scaled  *°£a 
of  the  walls  and  the  development  of  a farf leld  soiution  fo 
fixing  a boundary  condition  on  the  upstream  and  downstream 
boundaries  of  the  grid. 

The  treatment  of  upstream  and  downstream  boundaries 
follows  from  the  same  general  method  as  used  for  free  air  f o 
Ind  is  described  in  some  detail  in  Appendix  A.  As  summarized 
there,  an  asymptotic  solution  has  been  derived  by  the  method 
ofiSiqS  fTa  general  lifting  airfoil  in  a steady  subaonic 
flow  for  the  solid  wall  case.  This  solution  is  believed  to 
be  accurate,  and  is  also  used  for  the  case  of  small  wall 
porosities.  This  solution  is  thus  used  to  fix  a Dirichl 
houndarv  condition  for  <J>°  on  the  streamwise  boundaries  of 
the  finite  difference  grid  (i=l  and  i=IM) . No  corresponding 
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solution  for  the  unsteady  potential  has  been  developed;  how- 
ever the  steady  solution  (which  applies  for  the  unsteady  per- 
turbation  potential  for  k=0)  indicates  that  to  leading  order 
dx  is  exponentially  small  as  x-±».  Thus  a zero  gradient  con- 
dition d>  =0  is  imposed  on  i=l  and  IM  for  all  unsteady  cases. 

This  should  be  accurate  for  the  low  frequency  unsteady  flows 
considered  and  the  adequacy  of  the  approximation  has  been 
verified  for  the  quasi-steady  case.  Supersonic  wind  tu2n®?; 
flows  can  also  be  treated  by  the  method  and  in  this  case  the 
zero  gradient  upstream  and  downstream  condition  is  imposed 
for  both  steady  and  unsteady  solutions. 

Of  more  immediate  importance  to  the  numerical  simulation 
of  wall  interference  effects  is  the  treatment  of  the  wall  bound- 
ary condition  for  general  porous  or  slotted  walls.  The  solid 
wall  case  can  be  considered  as  a trivial  speciai  case  (Jy  0 
on  y = + h/2).  Analytical  modeling  of  ventilated  wind 
is  a subject  as  rich  in  problems  and  non-linear  effects  as  the 
transonic  flow  problem  itself.  Kacprzynski  has  in fa 
pointed  out  that  the  characteristics  of  porous  or  slotted  wal 
can  be  highly  dependent  on  flow  Mach  number  and  Reynolds  number 
in  addition  to  wall  construction.  Since  the  primary  purpose 
of  the  present  work  is  to  demonstrate  the  capability  of  the 
method  to  treat  wall  effects  it  was  decided  to  restrict  the 
initial  effort  to  the  simplest  possible  model.  a resu! It, 

the  wall  boundary  condition  is  treated  as  the  first  order  uin 
ear)  relation  between  the  perturbation  velocity  normal  to 

wall  and  the  local  pressure  coef £. icient,  as^^re  The  fSllow- 
m-essure  is  equal  to  the  freestream  static  pressure.  Tne  rono 

ing  equations  hold  for  the  steady  and  unsteady  perturbations 
respectively : 


4>: 


i ni  0 


= + P4> 


+ h 

°n  y “ - 2 


(42) 


<j>  = + PelB*  ( <|>  +ik  <t>) 


+ h 

on  y = - 2 


(43) 


where  the  scaled  tunnel  height  is. 


h = [(1+6)  «M<£]  /3  | 


(44) 


and  where  P is  the  wall  porosity  factor  (P  - 0 for  solid  walD 
and  3*  is  introduced  to  permit  a phase  lag  between  the 
pressure  and  normal  velocity  in  the  unsteady  case  The  bound- 
ary conditions  are  introduced  into  the  finite  ® JLr 

cedure  using  one  sided  difference  forms  for  4>yy  on  the  upper 


-41- 


and  lower  wall  boundaries  in  the  same  manner  as  the  body  bound- 
ary condition  and  as  described  for  the  supersonic  farfield 
condition  above. 

It  is  not  immediately  obvious  how  P is  to  be  estimated. 
Kacbrzvnski^^points  out  thit  it  could  be  a non-linear  function 
of  pressure  (or  »x>  and  that  P could  be  different  for  flow  into 
°r  ^ItTfthe  plenum.  In  any  event,  one  rational  estimate  for 

S23S  3 

ing  factors  are  the  open  area  ratio,  S,  or  1 S respective  y 
scaled  parameters: 


_ 6*^*  \J (j)  = 0 (free  jet)  (45) 

0,  1/3  X 

I (1  + Y>«*# 


and 


v = 6 U <J)y  = 0 


(solid  wall)  (46) 


where  u and  v are  the  streamwise  and  transverse  perturbation 
velocity  components.  Taking  the  weighted  average  gives = 


4>  = 


(1"s)  ki+y)sm^] 


TJ1  ^x 


(47) 


To  demonstrate  the  method,  steady  and  quasi:steady  (k-0) 
results  are  presented  in  Figures  l^an are  taken 

toCmatr?he“LR0PUot  ^n'nei»-since  the ;rfta  presented  in  Sec- 
tions 5.0  and  6.0  below  were  measured  there.  This  tunnel  n 
slotted  wall  with  a 10%  open  area  rating  tha^fo^th  ^ equa_ 

list: 

SiSgTf^e  ^cak.a£»”?rno?e%1haftheCp«rouTwfn9soluti^n 


-42- 


is  not  very  different  than  the  solid  wall  solution  indicating 
that,  for  this  case  at  least,  the  porosity  is  not  sufficient 
to  model  free  air  conditions.  Figure  18  presents  the  jump  in 
pressure  coefficient  per  unit  flap  angle  for  a quasi-steady 
(k=0)  flap  deflection  calculated  as  linearized  perturbations 
to  the  steady  solutions  just  described.  The  results  for  the 
solid  and  porous  wall  are  compared  to  the  corresponding  free 
air  solution  and  test  results  from  Reference  19.  The  results 
show  that  effect  of  the  walls  is  to  strengthen  the  perturbation 
through  the  shock  and  to  significantly  decrease  the  pressure 
perturbation  on  the  non-deflecting  forward  portion  of  the  airfoil. 
The  comparison  of  the  theoretical  results  show  that  the  wall  sig- 
nificantly decreases  the  upstream  influence  from  the  flap  to  the 
forward  portion  of  the  airfoil.  Although  the  wall  solutions 
do  not  improve  the  comparison  with  the  data  they  do  indicate 
that  appreciable  wall  effects  are  acting  for  this  test  config- 
uration. Additional  unsteady  results  with  wall  effect  are  pre- 
sented in  Section  6.0. 
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Figure  17.  Steady  Pressure  Coefficients 


O’ 


Three  Dimensional  Planar  Wing  Calculations 


M ,U 

oo  ^ 


Figure  19.  Schematic  of  Three-Dimensional  Planar  Wing 

The  application  of  the  theory  and  solution  method  to 
two-dimensional  airfoil  sections  presented  in  previous  work 
and  in  this  report  are  interesting  and  illustrative  but  for 
practical  application  to  dynamics  or  flutter  problems  three- 
dimensional  effects  must  be  considered.  As  with  most  °ther 
effects,  3-D  effects  are  more  important  at  transonic  speeds 
in  the  other  speed  ranges.  The  efficiency  of  the  present  scheme 
is  such  that  realistic  three-dimensional  computat ions  are  p a 
tical  on  modern  computers  and  it  is  the  Purpose  of  this  section 
to  describe  the  required  generalizations  and  to  demonstrate  the 
method  for  a simple  planform. 

The  wing  parameters  and  coordinate  definition  are  de- 
fined schematically  in  Figure  19  above.  The  initiai  develop- 
ment of  the  method  is  restricted  to  rectangular  Pljnformsun 
dergoing  oscillations  symmetric  with  respect  to  the  wing  root 
(z  = 0) . The  small  disturbance  analysis  and  the  unsteady  per 
turbation  theory  valid  for  three-dimensional  flows  were  de- 
scribed in  Section  2.  As  indicated  there,  the  generalization 
to  three  dimensions  requires  but  the  addition  of  the  term 

to  the  governing  equations  for  the  steady  and  unsteady  pertur 
bat ion  potentials.  Asymptotic  solutions  to  the  governing 
equations  have  been  derived  for  lifting  wings  in 
stream  flow  by  Klunker29,  for  the  steady  flow,  and  by  the  pre 
sent  authors  for  the  unsteady  perturbation.  These  solutions 
are  summarized  in  Appendix  A and  used  in  the  numerical  solution 
method  to  fix  farfield  boundary  conditions.  Three  dimensional 
solutions  for  steady  transonic  flow  have  been  presented  by 


Bailey  and  Steger30  and  Newman  and  Klunker'3'1';  the  latter  work 
being  most  closely  related  to  the  method  for  steady  flows 
used  in  this  work.  Extensions  of  the  solution  method  for  the 
unsteady  perturbation  parallel  the  steady  method  and  are  now 
described. 

3.4.1  Numerical  Solution  Method 


Figure  20.  Schematic  of  Numerical  Solution  Domain 

The  three-dimensional  numerical  scheme  constitutes  the 
most  straightforward  extension  of  the  two-dimensional  method 
previously  described  in  detail  in  Reference  1.  As  shown  sche 
matically  in  Figure  20,  a cubic  rectangular  mesh  of  finite  ex 
tent  with  uneven  grid  line  spacing  is  overlayed  on  the  3-D  so 
lution  space.  The  grid  is  concentrated  near  the  airfoil  and 
expanded  out  to  the  far  boundaries  of  the  grid.  The  finite 
difference  equations  are  identical  to  the  corresponding  two- 
dimensional  versions^  with  the  addition  of  a centered  differ- 
ence form  for  given  by: 
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Figure  21.  Schematic  of  Difference  Scheme 


The  computational  star  for  the  three  dimensional  scheme 
is  shown  schematically  in  Figure  21  above.  The  tests  for  the 
elliptic  on  hyperbolic  nature  of  the  equation  are  made  on  the 
centered  difference  form  of  ( K— ® ) < -i  an<^  depending  on  t e 
value  of  this  coefficient  the  x derivatives  of  <J>  are  centered 
or  backward  differenced  as  in  the  2-D  case.  A parabolic  point 
operator  is  used  in  both  the  steady  and  unsteady  schemes  and  a 
shock  point  operator,  as  described  for  the  2-D  case  in  Sec- 
tion  3.2,  is  used  in  the  finite  difference  scheme  for  the  stead* 

potential . 

As  before,  the  finite  difference  equations  are  set  up  for 
each  column  (x,z  = constants)  in  turn,  taking  account  of  the 


airfoil , wake  and  farfield 

nl^otTs 

The  linearization  is  accomplish  y ^ resuiting  linear  sys- 

ate  for  the  coefficient  V dx^y  optimuro  Gaussian  ellimin- 

tem  is  tri-diagonal  and  is  solv  Y £ terminated  when  the  dif- 
ation.  The  column  iterates  is  less  than  an  arbitrary 

ference  between  su^®S®i-5i  As  in  the  2-D  cane,  convergence 
small  amount  (usually  10  ).  As ^ iterations.  In  the  un- 

is usually  achieved  in  trailed  thftt  the  equation  is  linear  s 

that" noS co lumn” iteration  is  required. 

After  each  column  is  -^^dfo^thf tfcafnafure  of 
iable  relaxation  factor  which  ^Pf£oints  and  M „ .7  for  hyper- 
the  equation;  u>  ^ 1-7  for  el  nrncess  is  performed  for 

CtiCpoints.  The  left  to  right  in  x 

each  column  in  turn  sweeping  theg ^ farfield  (k  = KM)  in  z. 
and  from  the  wing  root  ( ^ edly  in  this  manner  until  the 

fhlngelr^for  a!lTia  pofnts  du^ing  one  grid  sweep  as  less 
than  some  arbitrary  small  amount. 

The  numerical  treatment  of teidfcasef if the^amf Is 

conditions  in  both thf  exception" that  the  airfoil  shape  fun c- 
the  2-D  case  with  Jhe  exceptio  x and  the  airfoil  cir 

tion  is  now  a function  of  * JJ  airfoil.  In  the  subsonic 

culation  is  a function  of  . £ asymptotic  solutions  for 

freestream  case  co£®ad®rgystems  described  in  Appendix  A are 
the  steady  and  unsteady  system  cQndition  on  five  sides  of 

used  to  fix  a D*ric^t  d Y containing  the  wing  root,.  a 
the  grid.  On  the  g 1 . s d whereby  <t2  = 0 on  z.  ' c 

?riUn"riXt”nn  depends  on  the 

SrffiSSS.'S  updatedSpef iodicallyCduring  the  solution  pro- 

C6SS • , i i 

The  solution  process  in- 

in  the  few  cases  calcula  e two-dimensional  met  o 

stance , seems  to  be  . »f  the  steady  and  un- 

steady^demonstration  calculations  Resented  in  the  next 

tion.  . - 

, Tnw  Asoect  Rectangular,  Zlau-g -- 

74?  Results  for  Low  Aspect  

J * ’ * ■ « • i-  i nrtC 
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To  demonstrate  the  three-dimensional  appUcsti°ns  o 

the  method,  results  are  presented  in  thrs 
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rectangular  planform  wing  with  a constant  six  percent  thick 
circular  arc  section  and  with  aspect  ratio  2b/c=8,  A grid 
of  approximately  12000  points  was  used  with  IM=30  over 
-3.2  < x < 3.4,  JM=19  over  -5.4  £ y < 5.4  and  KM=19  over 
0 < z~<  678.  In  the  x direction  16  grid  lines  were  distri- 
buted along  the  airfoil  with  Ax=.06,  and  in  the  z direction 
10  grid  lines  were  distributed  over  the  span  with  Az=.2. 

The  grid  design  was  chosen  so  that  the  unsteady  version  of 
the  program,  which  requires  three  times  the  grid  point  stor- 
age of  the  steady  program  (<j>°,  Re  <j>  and  Im  4>)  , would  fit  on 
small  core  of  the  CDC  7600  (^64K  storage). 

Results  for  the  steady  airfoil  pressure  distribution 
at  various  spanwise  stations  are  presented  in  Figure  22  be- 
low for  the  airfoil  described  above  at  a supercritical  Mach 
number  of  M^.86  and  zero  angle  of  attack.  The  correspond- 
ing two-dimensional  results  for  the  same  airfoil,  Mach  number 
and  (x,y)  grid  design,  are  shown  for  comparison.  It  is  noted 
that  the  pressure  distributions  for  the  inboard  50  percent  of 
the  semi-span  are  only  slightly  lower  than  the  two-dimensional 
solution  indicating  very  little  spanwise  relaxation  of  the 
flow  in  this  region.  The  pressure  does  not  begin  to  change 
significantly  until  the  final  25  percent  of  the  semi-span  near 
the  tip.  The  75  percent  station,  for  example,  shows  some  weak- 
ening and  forward  translation  of  the  shock  but  the  flow  does 
not  become  entirely  subcritical  until  very  near  the  tip  sta- 
tion. All  in  all,  the  solution  shows  a lesser  degree  of  three 
dimensionality  than  might  be  expected  for  this  case. 

Results  for  the  linearized  unsteady  perturbations  about 
the  steady  solution  just  described  are  presented  in  Figure  23 
and  24  for  a uniform  pitch  oscillation  at  k=0  and  k=  1,  re- 
spectively. In  each  case  the  corresponding  two-dimensional 
solution  is  also  shown.  These  results  for  the  jump  in  pressure 
across  the  airfoil  per  unit  pitch  angle  show  the  expected 
greater  degree  of  three  dimensional  effects  than  the  steady 
solution.  Of  particular  note  is  the  fact  that  the  root  sta- 
tion  in  each  case  is  considerably  different  in  magnitude  than 
the  2-D  solution.  This  is  especially  true  of  the  quasi-steady 
result  (Figure  23).  It  seems  that  the  main  effect  of  the 
three  dimensionality  is  in  depressing  the  peak  suction  through 
the  shock,  although  it  also  has  a considerable  effect  on  the 
phase  of  the  pressure  distribution  as  indicated  in  the  rela- 
tively  large  changes  in  Im(ACp~)  on  the  aft  portion  of  the  air- 
foil in  the  k=.l  case  (Figure  24).  It  should  be  noted  that 
the  results  show  some  residual  loading  on  the  airfoil  at  the 
tip  station  due  to  discretization  which  presumably  could  be 
reduced  by  using  a finer  mesh  near  the  tip.  Finally,  Figure  25 
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gives  the  span-wise  loading  distributions  for  each  of  the  3-D 
cases  with  comparison  to  the  corresponding  two-dimensional  re- 
sult. A comparison  of  the  Re(C£,a)  in  each  case  indicates  the 
important  result  that  the  reduced  frequency  effect  in  this 
low  frequency  regime  is  considerably  less  in  the  3-D  case  than 
in  the  2-D  case.  Such  a result  could  have  important  conse- 
quences for  a flutter  or  dynamic  stability  analysis  for  which 
the  variation  of  unsteady  force  coefficients  with  reduced  fre- 
quency can  have  a significant  effect  on  the  results. 

Some  comment  with  respect  to  computer  run  times  for 
the  cases  presented  here  is  clearly  of  interest.  The  steady 
computation  was  taken  to  a convergence  of  A<)>  °niax,v'5 . 0*10“ 
which  required  157  full  grid  iterations  and  was  performed  in 
a total  of  65  seconds  on  a CDC  7600.  The  quasi-steady  per- 
turbation required  286  grid  iterations  and  185  seconds  of  CDC 
7600  time  to  reach  a convergence  of  A<J>max=10~^.  Such  a fine 
degree  of  convergence  is  not  required  for  acceptable  accuracy 
and  it  is  noted  that  the  same  case  reached  A4>max=10“4  in  150 
iterations  and  82  seconds  of  7600  time.  Convergence  and  com- 
puter run  time  for  the  fully  unsteady  k=.l  case  were  similar 
to  (but  less  than)  the  k=0  case. 

Although  somewhat  greater  definition  (more  grid  points) 
is  called  for  than  was  used  in  these  demonstration  calculations, 
the  computer  time  requirements  mentioned  here  indicate  that 
the  method  is  efficient  enough  to  permit  flutter  aerodynamic 
calculations  for  three-dimensional  planar  wings.  The  three 
dimensional  versions  of  the  programs  which  implement  the 
method  discussed  above  for  the  steady  and  unsteady  perturba- 
tion potentials  are  described  in  some  detail  in  a companion 
users  manual  to  this  report.  The  method  as  presently  program- 
med applies  to  rectangular  planform  wings  in  subsonic  free- 
stream  flow  but  the  programs  could  be  easily  generalized  to 
treat  supersonic  or  wind  tunnel  flows.  The  generalization  to 
tapered  and  swept  wings  is  the  next  major  generalization  which 
could  be  effected  to  significantly  enhance  the  usefulness  of 
the  method. 
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Figure  23.  Jump  in  Airfoil  Pressure  Coefficient  Due  to 
Quasi-Steady  Deflection  in  Pitch 
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Figure  24.  Jump  in  Unsteady  Airfoil  Pressure  Coefficient 
Due  to  Oscillation  in  Pitch 


4.0  NUMERICAL  ACCURACY,  CONVERGENCE  AND  STABILITY 

CONSIDERATIONS 


Before  proceeding  with  the  large  number  of  calcula- 
tions presented  for  verification  of  the  method  (Section  5) 
and  for  the  flutter  aerodynamics  (Section  6)  the  numerical 
accuracy  and  convergence  of  the  method  was  investigated  in 
hopes  of  defining  a grid  and  degree  of  convergence  which 
would  strike  a happy  medium  between  accuracy  and  efficiency. 
Also,  attempts  were  made  to  speed  convergence  through  the 
optimization  of  relaxation  parameters  and  the  use  of  standard 
matrix  solution  acceleration  schemes.  These  studies  are  sum- 
marized in  Sections  4.1  and  4.2.  Studies  of  the  stability 
of  the  numerical  scheme  have  been  performed  throughout  the 
research  and  especially  upon  discovering  that  an  instability 
may  develop  for  a combination  of  freestream  Mach  numbers  near 
sonic  and  high  reduced  frequencies.  Some  salient  results  are 
summarized  in  Section  4.3. 


4. 1 Numerical  Accuracy  Considerations 


In  any  finite  difference  solution  technique,  numerical 
accuracy  is  often  as  important  as  the  accuracy  of  the  basic 
theory  upon  which  the  method  is  based.  The  ultimate  test  of 
numerical  accuracy  is,  of  course,  comparison  with  exact  solu- 
tions for  the  same  governing  equations.  Such  comparisons  have 
been  made  for  the  steady  solution  method15'16  and  these  have 
demonstrated  the  accuracy  of  the  mixed  differencing  line  re- 
laxation scheme  for  steady  flows.  Similar  comparisons  have 
been  made  for  the  present  unsteady  method  in  Section  4 for 
the  special  cases  of  linearized  subsonic  and  supersonic  flow 
and  have  indicated  that  the  solution  method  is  basically  ac- 
curate for  these  relatively  simple  flows.  It  is  unfortunate 
that  no  exact  solutions  are  available  for  unsteady  flow  at 
supercritical  Mach  numbers.  In  this  case,  however,  certain 
numerical  tests  can  be  performed  which  provide  a measure  of 
numerical  accuracy  and  some  results  of  such  tests  are  pre- 
sented here. 

The  accuracy  of  the  numerical  schemes  for  both  the 
steady  and  unsteady  perturbation  potentials  used  in  this  work 
depend  on  many  variables.  Two  of  the  most  important  are  the 
degree  of  convergence  of  the  iterative  relaxation  procedure 
and  grid  design  and  definition.  Their  effect  on  the  accuracy 
of  the  steady  and  unsteady  solution  methods  are  examined  in 
Figure  26  through  29  and  Figure  30  and  31  respectively. 
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Figure  26  shows  a typical  convergence  history  for  a super- 
critical steady  flow  calculation.  The  parameters  plotted 
versus  the  number  of  iterations  are  the  maximum  change  in  <J> 
throughout  the  grid,  (A4>max)  the  change  in  airfoil  circula- 
tion (Ayte)  and  the  change  in  the  non-uniform  flow  doublet 
distribution  (AP)  per  iteration,  where: 


V = 


(4>£) 2 d£dn 


(49) 


These  parameters  provide  a good  measure  of  the  degree  to  which 
the  solution  has  converged.  For  example  based  on  these  curves 
the  asymptotic  convergence  rate  for  each  parameter  can  be  de- 
fined and  used  to  estimate  the  error,  for  example  (<j>N-<J>oo), 
which  can  be  expected  if  the  iteration  procedure  is  terminated 
after  N iterations.  This  has  been  done  in  this  case  assuming 
the  iteration  procedure  was  terminated  when  A^a^lO-4  (N=*150) 
and  the  results  are  to  a conservative  estimate: 
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Thus  based  on  this  result  it  is  estimated  that  <j> 0 throughout 
the  grid  has  approached  its  limiting  value  to  within  0.11%; 
likewise  the  airfoil  lift.  This  is  clearly  more  than  accept- 
able accuracy  for  the  calculation  of  unsteady  aerodynamic 
coefficients  for  example;  the  prime  purpose  of  the  present 

work. 


Figure  27  presents  another  measure  of  accuracy  for 
this  same  case.  Presented  there  is  the  distribution  over  the 
airfoil  of  the  change  in  pressure  coefficient  due  to  change 
in  the  degree  of  convergence  from  A<J>max=10“  to  A4>max~10 
All  values  shown  in  the  figure  are  less  than  1%  of  the  cor- 
responding values  of  pressure  coefficient.  This  in  fact  is 
a more  vigorous  test  of  accuracy  than  considering  the  values 
of  <t>°  since  Cr/H^A^/Ax.  Thus  even  though  the  error  in  <f>  (A<|>) 
is  less  than  .1%  the  error  in  Cp  could  be  larger  due  to  the 
small  grid  size  on  the  airfoil.  This  figure  shows  the 
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interesting  fact  that  the  approximate  error  in  Cp  follows  the 
sairo  general  distribution  over  the  airfoil  as  the  pressure 
distribution;  with  the  error  being  greater,  as  miqht  be  ex- 
pected, near  the  leading  edge  and  at  the  shock. 

It  should  be  pointed  out  that  the  case  considered  is 
quite  typical  of  all  supercritical  cases  for  which  the  con- 
vergence has  been  examined  in  this  manner.  The  main  result 
of  such  investigations  is  that  a degree  of  convergence  of 
A^max^lO-^  should  be  more  than  sufficient  to  provide  1%  ac- 
curacy for  airfoil  pressures  and  integrated  forces. 

The  effect  of  grid  size  or  definition  on  solution  ac- 
curacy is  indicated  in  Figures  28  and  29.  Figure  28  presents 
results  for  the  pressure  distribution  for  a NACA64A006  airfoil 
at  M<»=. 875  using  two  different  grid  definitions:  a coarse  grid 

(25x25)  with  12  grid  lines  on  the  airfoil  and  a more  refined 
grid  (50x50)  with  26  grid  lines  on  the  airfoil.  Each  grid  was 
taken  to  the  same  degree  of  convergence  (A<j>max=10-4 ) and  the 
results  show  a relatively  large  difference.  The  comparison 
with  the  experimental  data  (Tijdeman  et  al!9,  see  Section  5) 
is  quite  good,  however,  for  the  refined  grid  so  that  no  fur- 
ther grid  refinement  was  used  in  this  case.  The  effect  of  an 
even  more  refined  grid  is  shown  in  Figure  29  for  a circular 
arc  airfoil.  As  indicated  there,  the  effect  of  such  a highly 
refined  grid  is  mainly  to  sharpen  up  the  shock.  The  effect 
on  the  pressure  distribution  away  from  the  shock  is  negligible. 
It  is  believed  that  the  same  sharpening  of  the  shock  could  be 
attained  by  concentrating  grid  lines  near  the  expected  shock 
position  rather  than  refining  the  entire  grid.  However,  in 
the  present  work  a grid  which  provided  acceptable  accuracy 
over  a wide  range  of  Mach  numbers  was  required  for  efficiency 
so  that  such  a strategy  for  sharpening  the  shock  could  not 
be  used  in  general.  Based  on  comparisons,  such  as  these,  it 
is  believed  that  a grid  with  approximately  50  points  in  each 
of  the  x and  y directions  will  provide  acceptable  accuracy. 


The  corresponding  effect  of  degree  of  convergence  on 
the  unsteady  scheme  is  shown  in  Figures  30  and  31.  The  case 
considered  is  the  NACA64A410  airfoil  at  M00=.  72  with  a quasi- 
steady pitch  deflection.  The  convergence  histories  for  A<j>max 
and  Ay  in  Figure  30  are  of  a somewhat  different  nature  than 
the  corresponding  steady  results  as  they  show  a relatively 
constant  rate  of  convergence  most  likely  due  to  the  linearity 
of  the  governing  equations.  Based  on  the  asymptotic  conver- 
gence history,  conservative  estimates  for  the  error  in  <(>  and  y 
if  the  iteration  process  is  terminated  at  A4>max=l°"’4  are: 
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Again  the  error  would  seem  to  be  acceptable  for  a degree  of 
convergence  of  A4>max=10“4  • The  corresponding  changes  m per- 
turbation pressures  on  the  airfoil  due  to  a cha^g®.^ rt°™ 
gence  from  A<t>max=10-4  to  A<)>max=10"5  are  shown  m Figure  31. 
The  distribution  is  shown  to  be  somewhat  erratic  but  in  e 
case  the  change  is  no  more  than  2%  of  the  local  value  of  Cp 
with  the  least  converged  pressures  being  at  the  leading 
edge  and  in  the  region  of  the  shock.  The  effects  of  gri 
definition  on  the  unsteady  perturbations  are  very  similar  to 
those  shown  for  the  steady  solution.  That  is,  the  major 
feet  of  a finer  grid  than  50  x 50  is  increased  definition  of 
perturbation  pressures  through  the  shock. 

In  summary,  studies  of  numerical  accuracy  as  briefly 
examined  here,  indicate  that  a 50  x 50  grid  taken  to  a con- 
veraence  level  of  A<j)max-10"4  provides  approximately  1%  ac 
curacy  in  both  steady  and  unsteady  results  for  pressures  on 
the  airfoil  and  a corresponding  level  of  numerical  accuracy 
for  integrated  forces. 
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Figure  31.  Influence  of  Degree  of  Convergence  on 
Unsteady  Airfoil  Pressures 
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Convergence  of  the  steady  and  unsteady  solution  schemes 
is  influenced  by  the  parameters  which  define  the  solution 
(similarity  parameters  K,ft,  airfoil  shape  function)  and  by 
the  parameters  which  control  the  numerical  solution  procedure 
(grid  design,  and  relaxation  parameters).  Based  on  the  large 
number  of  solutions  calculated  for  this  report,  using  the 
basic  50  x 50  grid,  the  convergence  histories  presented  above 
are  typical  for  a wide  range  of  supercritical  Mach  numbers  and 
reduced  frequencies.  The  relaxation  parameters  can  have  a 
somewhat  greater  effect  on  convergence  and  this  effect  is  dis- 
cussed here  along  with  a discussion  of  possible  convergence 
acceleration  schemas. 

As  noted  in  previous  work,  a type  dependent  relaxation 
scheme  is  used  in  both  the  steady  and  unsteady  solutions;  that 
is,  for  locally  subsonic  points  the  solution  is  over-relaxed, 
l<o)e<2  and  for  locally  supersonic  points  the  solution  is  under- 
relaxed, 0<o)h<i.  The  limits  on  the  relaxation  parameters  are 
due  to  stability  considerations  and  are  common  for  elliptic 
or  hyperbolic  finite  difference  solution  schemes.  The  values 
of  the  relaxation  parameter  within  the  above  limits  can  affect 
the  convergence  rate  and  in  general  there  exist  optimum  values 
which  provide  the  fastest  convergence.  The  effect  of  the  el- 
liptic relaxation  factor,  we,  on  convergence  for  a typical 
supercritical  steady  case  is  shown  in  Figure  32.  It  shows 
that  for  a value  of  u»e=1.7,  A(J>max=10_4  convergence  is  reached 
in  58  iterations  as  compared  to  78  for  o)e=1.5.  Almost  identi- 
cal results  as  shown  for  we=1.7  were  obtained  with  o)e=l»8  and 
for  a value  of  o)e=1.9  the  convergence  was  somewhat  faster  but 
was  oscillitory  in  nature.  Such  oscillations  would  indicate 
that  a)<v,1.7  or  1.8  is  optimum.  The  effect  of  the  hyperbolic 
relaxation  parameter  in  this  case  and  for  all  subsonic  cases 
considered  was  negligible  for  w ^=.7  to  .9.  A similar  result 
holds  for  the  unsteady  solutions  examined  parametrically  in 
this  same  fashion.  Thus  values  of  to  1.8  and 

to  .9  seem  to  be  close  to  optimum  for  a wide  range  of  simi- 
larity parameters  of  interest.  In  other  words,  for  any  given 
case  the  convergence  rate  could  probably  be  increased  by 
modest  amounts  but  the  experimentation  needed  to  determine 
optimum  values  would  outweigh  the  increased  efficiency  achieved. 
The  task  of  mapping  out  optimum  values  of  we  and  u>h  as  func- 
tions of  the  similarity  parameters  K and  ft  would  be  an  under- 
taking which  would  not  be  justified  due  to  the  surprisingly 
small  effect  the  parameters  seem  to  have  for  values  in  the 
range  given  above.  This  is  especially  true  since  any  bene- 
fits so  achieved  might  possibly  be  negated  by  grid  geometery 
or  airfoil  shape  effects. 


For  supersonic  freestream  flows  the  importance  of  we 
relative  to  oih  switches.  That  is,  becomes  the  more  im- 
portant relaxation  factor.  However,  based  on  the  supersonic 
calculations  performed  in  this  study  the  ranges  of  these 
parameters  given  above  still  provide  close  to  optimum  results. 
It  should  be  mentioned  that  for  fully  supersonic  flow  the 
present  scheme  converges  in  one  grid  iteration,  since  in 
this  case  the  scheme  provides  an  implicit  solution  method  for 
the  wave  equation  in  which  x is  the  time  like  variable.  This 
is  true,  however,  only  if  w^=l.  For  w^l  in  the  fully  super- 
sonic case,  the  scheme  may  require  many  tens  of  iterations 
to  converge. 

The  possibility  of  increasing  the  efficiency  of  the 
solution  mehtod  by  using  various  matrix  solution  accelera- 
tion schemes  has  been  investigated  briefly  in  this  study. 

As  described  in  previous  work,1'20  one  such  method  which  has 
worked  well  is  the  so-called  "grid  halving"  procedure  in 
which  the  solution  is  initiated  on  a coarse  grid,  taken  to 
a certain  level  of  convergence,  and  successively  refined. 

Another  method  which  was  tried  was  the  Aitken's  6 2 
method33.  This  method  in  its  most  basic  form  assumes  that 
a sequence  of  numbers  (4>n”  ,<|>n,(J>n  / converges  exponentially 

with  respect  to  n and  is  used  successfully  to  accelerate 
convergence  of  iterative  schemes  for  system  of  equations. 

For  use  in  the  present  scheme,  on  extrapolated  value  for  the 
potential  at  any  field  point  can  be  defined  based  on  three 
successive  iterates  by: 
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(50) 


In  the  present  work  this  has  been  applied  in  two  ways:  the 

first  is  to  extrapolate  all  grid  points  and  the  second  is 
to  extrapolate  a global  quantity  such  as  airfoil  circulation. 

The  use  of  extrapolation  on  all  grid  points  for  the 
non-linear  steady  potential  has  given  quite  bad  results  in 
all  cases  tried.  The  non-linearity  of  the  basic  equations 
is  a possible  culprit  but  similarly  unpleasant  results  were 
obtained  when  the  method  was  applied  to  the  linear  unsteady 
perturbation  potential.  The  use  of  extrapolation  on  the  air- 
foil circulation  has  provided  mixed  results;  that  is  it  has 
worked  well  in  some  cases  but  poorly  in  others.  However, 
even  in  the  cases  in  which  it  worked,  the  overall  convergence 
of  the  scheme  was  little  if  any  better  than  the  use  of  over-  34 
relaxation.  A higher  order  extrapolation  scheme  after  Shanks 
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which  generalizes  the  Aitken  method  was  also  tried  with 
similar  bad  results. 

As  a postscript , it  should  be  mentioned  that  the  sing- 
ular lack  of  success  with  acceleration  methods  achieved  during 
this  study  should  not  discourage  further  investigations  along 
these  lines  because  of  the  possible  significant  increases  in 
computational  efficiency  which  such  schemes  promise.  In  tact, 
Hafez  and  Cheng^  have  recently  demonstrated  that  accelera- 
tion can  be  made  to  work  even  for  the  non-linear  transonic 
potential  equation.  Their  method,  however,  was  not  available 
in  time  to  incorporate  into  the  present  work. 
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4 . 3 Stability  Considerations 


Stable,  convergent  steady  and  unsteady  calculations 
have  been  performed  for  a wide  range  of  Mach  numbers  and 
reduced  frequencies  as  testified  to  by  the  numerous  results 
presented  in  this  report.  Stability  problems  have  however 
been  encountered  with  the  steady  and  unsteady  solution  schemes 
for  near  sonic  freestream  Mach  numbers. 

In  the  present  work  instabilities  have  developed  for 
both  steady  subuonic  and  supersonic  freestream  calculations 
near  Mach  1.  Such  problems  have  been  alluded  to  by  previous 
investigators^  and  Kentzer-^^  has  shown  the  potential  for  a 
long  wave  instability  in  the  Murman-Cole  scheme.  The  insta- 
bility does  not  seem  to  be  of  practical  importance,  however, 
since  many  authors  have  applied  the  same  basic  scheme  to 
steady  flows  with  Mach  numbers  quite  close  to  one.  With  re- 
spect to  the  problems  encountered  in  the  present  work  it 
would  seem  that  they  are  more  likely  due  to  inadequate  grid 
design.  For  example  the  scheme  becomes  unstable  if  super- 
sonice  flow  regions  approach  the  grid  boundaries  for  subsonic 
freemstream  flow  cases  and  if  regions  of  subsonic  flow  approach 
the  grid  boundaries  for  supersonic  freestream  flow.  Such 
problems  can  usually  be  solved  by  expanding  the  grid  farther 
from  the  airfoil.  In  the  calculations  reported  here  the  basic 
grid  with  -6<y<6  could  be  used  for  the  NACA64A006  up  to  M^.95 
and  the  NACA64A410  to  M =.9.  In  the  supersonic  flow  calcula- 
tions for  the  NACA64A410  it  was  necessary  to  expand  the  grid 
to  -12<y<12  to  perform  calculations  for  freestream  Mach  num- 
ber as  low  as  Mw=l.l.  A comprehensive  study  of  these  problems 
was  beyond  the  scope  of  the  present  work. 

The  solution  method  for  the  unsteady  potential  was 
found  to  have  a more  restrictive  stability  problem  for  near 
sonic  flows.  Stability  problems  have  been  encountered  for 
a combination  of  Mach  numbers  near  1 and/or  increasing  re- 
duced frequency.  For  example  for  the  NACA64A006  airfoil  at 
M^.95  the  method  was  unstable  for  k> . 1 and  for  all  reduced 
frequencies  at  M =1.05. 

In  previous  work  it  was  noted  that  the  scheme  could  have 
stability  problems  near  a compression  or  shock.  A heuristic 
stability  analysis  indicates  that  the  instability  is  related  to 
the  inability  of  the  finite  difference  grid  to  resolve  short 
wavelength  (X„)  upstream  facing  waves.  That  is,  the  requirement 
thc.c  Au=  (1-M)  /k>Ax  is  not  satisfied  as  M (local  Mach  number) 
approaches  1 and  k increases.  If  this  occurs  over  a large  enough 
region  of  the  grid  an  instability  is  likely.  It  is  noted  that 
expansion  of  the  grid  as  used  here  most  likely  compounds  the 
problem. 


\ 
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A more  rigorous  Von  Neumann  stability  analysis  has 
also  been  performed  for  the  present  finite  difference  scheme. 
The  resulting  stability  expression  is  quite  complicated  de- 
pending on  Ax,  Ay,  k,  ft,  we,  4>x'  4>xx  as  well  as  the  usual 

disturbance  wavelength,  and  it  need  not  be  given  here.  For 
the  simple  case  of  uniform  flow  (<j>x=4>xx=0)  ' the  expression 
indicates  that  the  scheme  is  stable  for  all  reduced  fre- 
quencies with  the  usual  restrictions  on  the  relaxation  fac- 
tors. For  general  supercritical  flows,  the  stability  ex- 
pression indicates  the  potential  for  instability  for  small 
values  of  K(M00->1)  for  large  values  of  ft  (or  k)  and  for  long 
wavelength  disturbances.  This  is  similar  to  Kentzer’s  re- 
sult for  the  steady  scheme  and  it  means  that  disturbances 
(due  to  round  off  etc,)  with  a wavelength  long  in  relation 
to  the  grid  size  can  amplify  and  therefore  destroy  the  solu- 
tion. A similar  description  of  this  instability  has  been 
developed  by  Ehlers-^  using  an  eigenvalue  analysis  of  the 
column  relaxation  scheme.  It  is  noted  that  this  instability 
works  at  cross  purposes  with  the  requirement  for  a grid  of 
ever  increasing  extent  as  flow  Mach  number  approaches  sonic. 
Increasing  the  size  of  the  grid  would  in  fact  enhance  the 
likelihood  that  this  long  wavelength  instability  would  develop. 

The  problem  can  not  probably  be  resolved  by  a simple 
change  to  the  finite  difference  equations.  One  possible  sol- 
ution would  be  to  use  the  results  of  an  inner/outer  expansion 
procedure  which  suggests  that  near  the  airfoil  the  present 
formulation  (Equation  10)  applies  but  away  from  the  airfoil 
the  high  frequency  equation  (Equation  21)  applies  if  k>0(M-l) 
where  M is  the  local  flow  Mach  number.  Thus  near  the  airfoil 
where  the  local  Mach  number  may  be  significantly  different 
from  sonic  the  present  finite  difference  procedure  could  be 
used  and  matched  to  an  available  analytical  solution  for  the 
high  frequency  equation  away  from  the  airfoil  where  local 
Mach  number  approaches  freestream.  The  regimes  of  flow  Mach 
number  and  reduced  frequency  for  which  this  procedure  is 
valid  would  be  somewhat  restricted  but  it  would  most  surely 
reduce  the  size  of  the  required  finite  difference  grid  for 
M^l  thereby  potentially  elliminating  the  stability  problems. 


5.0  COMPARISONS  WITH  ALTERNATE  THEORIES  AND  EXPERIMENTS 


The  accuracy  of  small  disturbance  theory  for  %-eady 
inviscid  transonic  flows  has  been  well  documented  in  recent 
years  through  comparison  with  exact  solutions,  detailed  nu- 
merical calculations  and  experiments  in  which  viscous  effects 
are  small.  The  purpose  of  this  section  is  to  provide  some 
initial  insight  into  the  potential  accuracy  of  the  present 
small  perturbation  method  for  unsteady  transonic  flows.  In 
previous  work-1- it  has  been  shown  that  the  present  approach 
matches  linearized  unsteady  subsonic  theory  for  subcritical 
Mach  numbers  and  the  results  of  Section  3.2  above  show  that 
the  method  matches  linearized  supersonic  theory  for  fully 
supersonic  flows.  These  trends  are  verified  by  the  aerody- 
namic coefficients  presented  in  this  section  and  the  next. 
Although  no  exact  solutions  exist  for  unsteady  transonic  flows, 
Magnus  and  Yoshihara^  have  recently  presented  a detailed  nu- 
merical calculation  and  it  is  compared  to  the  present  theory 
in  Section  5.1.  The  final  verification  of  any  theory  must  of 
course  come  from  comparison  to  experimental  data  and  such  a 
comparison  is  made  in  Section  5.2. 

Results  are  presented  for  a NACA  64A410  airfoil  oscil- 
lating in  pitch  and  a symmetric  NACA  64A006  with  oscillating 
control  surface  hinged  at  the  3/4  chord  line.  Before  proceed- 
ing with  the  comparisons,  some  comments  concerning  computa- 
tional details  are  in  order. 

In  all  of  the  unsteady  calculations  to  follow,  the  low 
frequency  approximation,  which  involves  dropping  underlined 
terms  in  the  above  formulation,  was  used.  The  effect  that 
this  approximation  might  have  on  the  results  has  been  discussed 
in  Section  3.1  above.  The  form  of  the  unsteady  perturbation 
to  the  airfoil  downwash  for  pitch,  plunge,  and  control  surface 
oscillations  are  respectively: 
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where  xh  is  the  hinge  point  and  xD  is  the  pitch  point.  In  the 
respective  cases  , the  magnitude  of  the  perturbation  is  ot/o  , ( / 
or  6 ' /6  where  a,P>  and  6'  are  the  amplitudes  of  the  oscillating 
component  of  angle  of  attack,  flap  angle  or  plunge  displacement 
respectively.  The  amplitudes  of  the  resulting  perturbations  to 
the  airfoil  pressure  coefficients  in  unsealed  (Cp)  or  scaled  (Cp) 
form  are  given  by  Equations  11  and  12  above.  It  is  noted  that 
the  scaling  used  in  Equation  11  for  the  pressure  coefficient, 
commonly  known  as  Spreiter  scaling,  is  a natural  result  of  the 
non-dimensionalization  used  here.  No  attempt  has  been  made  to 
modify  the  Mach  number  dependence  in  the  scaling  to  provide 
better  aareement  with  data  or  exact  results  as  is  often  done 
elsewhere16.  Finally,  it  is  noted  that  lift,  moment,  and  hinge 
moment  coefficients  (Ex.  Cect,  Cmct,  Cna)  are  calculated  by  a trap- 
eziodal  rule  integration  of  the  calculated  pressure  coefficient 
distributions . 


In  the  calculations  presented  in  this  section  and  the 
next,  the  refined  grid  (achieved  by  refining  a coarse  grid  by 
the  "grid  halving"  procedure)  consisted  of  48  points  the 
x direction  and  47  points  in  the  y direction.  Over  half  (26) 
of  the  grid  columns  are  on  the  airfoil  and  grid  rows  are  di- 
vided in  a symmetric  manner  above  and  below  the  airfoil.  The 
qrid  is  expanded  in  a regular  manner  out  to  x = ±3  and  y - 
at  which  point  an  analytic  farfield  expression  is  ^p°sed  as 
a boundary  condition.  Previous  experience  indicated  tha*  auch 
a grid  would  provide  relatively  good  results  so  that  no  attempt 
was  made  at  optimizing  the  grid  structure  for  each  case, 
both  steady  and  unsteady  calculations  the  iteration  procedure 
was  continued  until  the  maximum  change  in  potential,  over  ail 
grid  points,  from  one  iteration  to  the  next  was  less  than 
X 10”*.  This  is  consistent  with  the  accuracy  tests  described 
in  Section  4.0  above.  In  each  case  this  resulted  in  a conver- 
gence of  the  airfoil  circulation  to  better  than  10  . 


5.1  NACA  64A410  Oscillating  in  Pitch  at  Ma 


= .72 


Magnus  and  Yoshihara9  (M-Y)  have  recently  presented 
time-dependent  calculations  of  the  flow  past  a NACA  64A410 
airfoil  oscillating  in  pitch  at  Mach  .72.  Their  highly  re- 
fined numerical  method  is  based  on  the  full  unsteady  inviscid 
flow  equations  with  kinematic  boundary  conditions  fixed  at  the 
mean  airfoil  surface.  These  "exact"  solutions  provide  a unique 
comparison  for  verification  of  the  present  approximate  method. 
ThePnose  bluntness,  thickness  and  camber  of  the  J;?upled 

with  the  mean  and  oscillating  amplitude  of  angle  of  a£tack 
(a  “ 2°  ± 2°)  provide  a severe  but  practical  test  of  botn  tne 
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small  disturbance  approximation  and  the  linearized  unsteady 
method. 

Results  for  the  steady  pressure  distributions  for 
a0  = 2°  are  presented  in  Figure  33  below;  pluses  (+)  and  A 
indicate  points  computed  by  the  present  method  and  M-Y  res- 
pectively. All  in  all  the  comparison  is  moderately  good  with 
the  exception  of  two  anomalies;  the  shock  description  and  a 
"knee"  in  the  small  disturbance  pressure  distribution  on  the 
upper  surface  near  the  nose.  Lift  predicted  by  both  methods 
compare  well  as  do  the  lower  surface  pressure  distributions. 

The  shock  calculated  by  the  present  method  is  farther  forward, 
more  diffuse  and  weaker  which  can  be  attributed  to  the  rela- 
tively coarse  grid  and  the  use  of  the  non-conservative  differ- 
encing of  the  present  method.  It  is  believed  that  the  shock 
description  could  be  improved  by  the  use  of  shock  point  dif- 
ferencing^ described  in  Section  3.2  or  shock  fitting^  as  well 
as  a denser  gridwork  in  the  region  of  the  shock.  However,  it 
should  be  pointed  out  that  the  shock  Mach  number  is  about  1.3 
which  would  indicate  that  the  small  disturbance  shock  descrip- 
tion is  not  accurate.  The  other  anomaly  in  the  comparison 
which  has  some  bearing  on  the  linear  perturbation  results  to 
follow  is  the  sharp  change  in  slope  of  the  predicted  pressure 
distribution  near  the  nose.  This  behavior  is  not  shown  in 
the  M-Y  results.  It  is  probably  due  to  the  approximation  to 
the  airfoil  slope  distribution  used  in  the  present  method  as 
well  as  possible  inaccuracy  of  the  small  disturbance  descrip- 
tion of  the  rapid  expansion  of  the  flow  about  the  nose  of  the 
airfoil.  It  is  clear  from  this  comparison  that  the  present 
results  could  benefit  from  increased  resolution  as  well  as 
the  use  of  recent  developments  in  the  treatment  of  shocks. 

Calculations  of  linearized  unsteady  perturbations  about 
the  steady  flow  described  above  have  been  performed  and  ampli- 
tudes of  perturbations  to  the  pressure  coefficents  due  to  a 
2°  pitch  oscillation  are  shown  in  Figures  34  and  35  for  the 
steady,  k=0,  and  low  frequency  unsteady,  k=.2,  cases.  The 
results  of  the  present  theory  are  again  compared  to  results 
deduced  from  the  M-Y  calculations.  The  M-Y  results  for  both 
nose  up  (2°-*'4°)  and  nose  down  (2°+0°)  motion  are  shown  to  give 
some  indication  of  the  inherent  nonlinearity.  It  is  noted  that 
the  nonlinear  behavior  is  surprisingly  small  on  the  super- 
critical upper  surface  and  surprisingly  large  on  the  sub- 
critical  lower  surface.  The  present  results  for  the  lower 
surface  compare  well  with  the  M-Y  results.  The  comparison 
in  each  case  for  the  upper  surface  is  consistent  with  the 
difficiencies  in  the  present  result  for  the  steady  flow  being 
perturbed.  The  discrepancy  at  the  shock  is  to  be  expected 
since  the  M-Y  calculations  show  significant  shock  excursions 
(indicated  on  the  graphs)  which  are  not  accounted  for  in  the 


present  work.  The  linearized  results  show  a peak  suction 
through  the  steady  shock  which  is  an  attempt  by  the  linear 
theory  to  account  for  the  strengthening  (or  weakening)  of  the 
shock  due  to  a nose  up  (or  nose  down)  pitch  perturbation.  The 
present  results  also  show  an  as  yet  unexplained  decreased 
suction  (dip)  near  the  x=.l  station  which  corresponds  to  the 
"knee"  in  the  steady  pressure  distribution  discussed  above. 
This  behavior  is  not  shown  in  the  M— Y results  and  is  an  indi- 
cation of  possible  numerical  difficulties  in  the  linear  per- 
turbation solution.  This  is  not  believed  to  be  the  case, 
however,  since  a comparison  of  nonlinear  steady  small  distur- 
bance results  (not  shown)  for  oo=0°  and  <Xq=2  show  a similar 
dip  in  the  distribution  of  ACpfl->2°.  This  would  indicate  that 
the  anomaly  is  due  to  the  discrepancy  in  the  steady  solution 
in  this  same  region  as  discussed  above. 
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Table  1.  Comparison  of  Lift  and  Moment  Pertrubations 


The  results  of  interest  for  flutter  applications  are 
the  perturbations  to  the  integrated  airfoil  forces.  The  table 
presented  above  compares  the  increment  in  lift  and  moment  coef- 
ficient in  steady  flow  as  the  angle  of  attack  increases  from 
0°  to  2°  and  from  2°  to  4°.  Also  compared  are  the  linear  parts 
of  the  oscillating  coefficients  for  the  k=.2  pitching  oscillation 
with  2°  amplitude.  The  good  agreement  in  each  case  would  seem 
to  indicate  that  the  linearized  perturbation  method  provides 
a useful  approximation  to  the  oscillatory  forces  for  quasi-steady 
and  low  frequency  oscillations.  It  may  be  noted  that  for  this 
example  the  ratio  of  oscillatory  amplitude  (at  2°)  to  the  thick- 
ness is  e=.  35  which  is  a relatively  large  perturbation. 


I 


I 

The  M-Y  results  for  k=.2  shown  in  the  table  are  the 
magnitudes  of  the  first  harmonic  (eiwt)  coefficients  which 
resulted  from  their  Fourier  analysis  of  the  nonlinear  time 
dependent  results.  Their  analysis  also  shows  that  the  higher 
harmonics  are  quite  small  in  comparison.  This  is  to  be  ex- 
pected from  the  perturbation  expansion  which  shows  that  terms 
involving  higher  harmonics  are  of  order  e2.  This  result  lends 
some  credence  to  the  linearized  approach  used  here.  It  is 
hoped  that  future  improvements  in  both  the  small  disturbance 
steady  and  linearized  unsteady  perturbation  methods  will  im- 
prove certain  details  of  the  present  results. 
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Figure  34.  Magnitude  of  Pressure  Coefficient  Perturbation 
Due  to  a 2°  Deflection  in  Pitch 
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Magnus  & Yoshihara 
Nose  Up  a=2°->4° 
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Figure  35.  Magnitude  of  Pressure  Coefficient  Pertur 
bat ions  Due  to  a 2°  Oscillation  in  Pitch 
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NACA 


64A006  with  Control  Surface  Oscillation 


As  mentioned  in  the  introduction,  the  primary  applica- 
tion of  the  method  described  here  is  the  calculation  of  gener 
alized  unsteady  forces  for  use  in  flutter  studies.  Such  a 
study  has  been  undertaken  for  the  NACA  64A006  airfoil  wit 
three  rigid  body  degrees  of  freedom:  pitch,  plunge  and  control 

surface  oscillation;  and  the  detailed  results  are  presented 
in  Section  6.0.  In  this  section,  some  of  the  results  generated 
for  control  surface  oscillations  are  Presented  with  comparison 
to  the  experimental  data  of  Tijdeman,  Bergh  ^ *ppe 

for  the  same  airfoil  and  control  surface  configuration.  It 
is  noted  that  the  experimental  technique  used  to  measure  un- 
steady pressures  was  sensitive  only  to  the  first  har")°"1^ 
so  that  theory  and  experiment  are  being  compared  on  the 

bases . 

Sample  results  for  airfoil  pressure  coefficients  (in 
scaled  form)  for  the  non-linear  steady  (Figure  36)  and  for 
linearized  perturbations  due  to  control  surface 
for  k=0  (Figure  37)  and  k=.12  (Figure  38)  are  presented  below. 

The  steady  results  for  a supercritical  Mach  number  °f 

and  zero  degree  angle  of  attack  show  good  agreement  with  eX 

Deriment  with  the  exception  of  the  region  near  the  shock. 

The  present  small  disturbance  results  show  what  must  be  called 
a rapid  compression  at  about  midchord  whereas  the  data  shows 
a diffuse  but  moderately  strong  shock.  Refinement  of  the 
grid  between  x/c=.5  and  .6  would  most  surely  enhance  the  com 
parison  but  it  is  likely  that  the  use  of  point  dif ferenc 

ing  would  over  correct  for  the  rearward  position  of  the  shock. 

The  linearized  perturbation  results  for  a steady  con- 
trol surface  deflection  and  for  an  unsteady  control  surface 

oscillation  at  k=.12  (Figures  37  and  38 ) show tj^theo™  over- 
agreement  with  data.  As  shown  in  Figure  37,  the  theory  over 

predicts  the  pressure  perturbation  on. th®,  st®^°narL  • 
nortion  of  the  airfoil  in  comparison  to  the  data.  Tnis  is 
characteristic  of  all  cases  calculated  and  is  believed  to  be 
due  in  part  to  the  reduced  upstream  influence  *e Trls- 
wall  effects  in  the  experiment.  It  is  ^a^ledc^a^lati0nS 
ent  calculations  are  for  free  air  conditions . 
with  a solid  wall  boundary  condition,  presented  in  Section 
614,  show  a significant  reduction  in  upstream  influence. 

Other  details  of  the  pressure  distribution  such  as  peak  suc- 
tion at  the  shock  (slightly  forward  of  the  by the 

hinae  Doint  singularity  seem  to  be  well  predicted  by  tne 

theory  9 The  unsteady  k=.12  results  in  Figure  28  show  a similar 

qualitative  comparison  with  some  significant  qualitative  dif 
ferences  in  magnitude  and  phase  on  the  forward  portion  of  the 

airfoil . 
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Comparisons  of  the  theoretical  results  for  with 
experimental  data  and  with  linear  subosnic/supersomc  theory 
are  made  in  Figures  39  and  40  for  k=0  and  0.1  respectively. 

As  shown,  the  present  theory  reduces  to  linear  theory  away 
from  Mach  1.  In  each  case  the  present  results  are  signifi- 
cantly greater  than  the  data  although  a qualitative  agreement 
is  apparent.  The  peak  in  the  theoretical  results  at  M=.85 
seems  to  result  from  a lift  augmentation  effect  due  to  coal- 
escence at  the  shock  of  upstream  traveling  control  surface 
perturbations.  This  effect  does  not  occur  in  linear,  uniform 
flow  theory  for  obvious  reasons  and  is  not  in  evidence  in  the 
data  perhaps  because  of  the  reduced  upstream  influence  men- 
tioned above.  One  final  comment  with  respect  to  these  compar- 
isons is  in  order.  It  is  noted  that  the  data  is  in  considerable 
disagreement  with  subsonic  linear  theory  (and  the  present 
theory)  for  a fully  subcritical  Mach  number  of  .7.  Nonuniform 
flow  or  wind  tunnel  wall  effects  could  not  be  the  culprits 
because  they  both  have  an  effect  which  would  tend  to  enhance 
the  discrepancy.  It  seems  highly  probable  that  the  reasons  for 
the  discrepancy  is  the  existence  of  a strong  viscous  effect 
in  the  data.  The  estimated  chord  Reynold's  number  of  less  than 
200,000  for  the  experiments  adds  further  weight  to  this  con- 
clusion. Such  effects  could  only  be  magnified  in  the  transonic 
speed  range. 


Further  comparisons  with  the  data  are  presented  in  the 
next  section  along  with  some  results  for  mean  angle  of  attack 
and  wall  effects. 
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Figure  39.  Lift  Due  to  Aileron  Deflection 
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Figure  40,  Lift  Due  to  Aileron  Oscillation 


6.0 


CALCULATED  RESULTS 


Using  the  method  described  in  this  report,  a wide 
range  of  steady  and  unsteady  results  have  been  calculated 
to  provide  the  necessary  aerodynamic  coefficient  data  for 
the  flutter  studies  to  be  presented  in  Section  7.  These  re- 
sults, described  in  this  section,  are  for  two  practical  air- 
foil: a NACA64A006  oscillating  in  pitch,  plunge  and  control 
surface  rotation  and  a NACA64A410  oscillating  in  pitch  and 
plunge.  The  ranges  of  Mach  number  (Moo=.75  to  1.2)  and  re- 
duced  frequency  (k=0  to  0.2)  considered  serve  to  demonstrate 
the  efficiency  and  utility  of  the  method  for  treating  un- 
steady transonic  flows. 


The  same  basic  grid  as  just  described  in  Section  5 
(IM=JM-50,-3<x<3, -6<y<6)  was  used  in  all  cases  with  but 
slight  variations  for  some  special  cases.  Also,  as  is  the 
usual  practice  with  the  method,  convergence  in  each  case  was 
taken  to  }.,-jA<f>*10  4,  which  resulted  in  the  convergence  of 
airfoil  circulation  to  better  that  10“4.  it  is  reiterated 
that  this  grid  design  and  degree  of  convergence  represents 
a reasonable  compromise  between  numerical  accuracy  and  com- 
puter run  time. 


6*1  Results  for  a NACA  64A006  Airfoil 


Results  are  presented  in  this  section  for  a NACA64A006 
airfoil  with  a 25%  trailing  edge  flap.  Detailed  steady  and 
unsteady  pressure  distributions  are  presented  in  Section  6.1.1 
and  the  calculated  unsteady  aerodynamic  coefficients  for  pitch, 
plunge  and  control  surface  oscillations  are  summarized  in  Sec- 
tion 6.1.2.  Mean  angle  of  attack  and  wall  interference  effects 
on  the  aerodynamic  coefficients  are  briefly  examined  in  Section 
6.1.3  and  6.1.4  respectively.  It  is  noted  that  the  unsteady 
results  for  pitch  and  control  surface  oscillations  were  calcu- 
lated using  the  low  frequency  approximation  since  the  high  fre- 
quency generalizations  to  the  method,  required  and  used  for 
the  plunge  results,  were  not  inplemented  at  the  time  these  cal- 
culations were  performed.  The  effect  of  the  approximation  on 
the  results  has  been  discussed  in  Section  3.1. 


6*1*1  Free  Air  Results  of  Steady  and  Unsteady  Pressures 


Steady  pressure  distributions  for  the  NACA64A006  at 
zero  angle  of  attack  for  subsonic  freestream  Mach  numbers, 
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Moo®. 8,  .85,  .9,  .95,  and  for  supersonic  freestream  Mach  num- 
bers, Mob-1. 1,  1.15,  1.2,  are  presented  in  Figure  41  and  42 
respectively.  The  various  curves  show  the  progression  of 
airfoil  pressures  from  subcritical  through  the  supercritical 
regime  to  low  supersonic  and  points  out  the  large  non-linear 
effect  in  the  transonic  regime.  The  airfoil  first  becomes 
supercritical  at  a Mach  number  between  .8  and  .85  and  by  M^— 

.85  a weak  shock  has  formed  which  strengthens  and  moves  aft 
on  the  airfoil  as  M<»  increases.  At  Mco=.9  the  shock  is  about 
at  the  aileron  hinge  point  (x/c=.75)  and  at  Moo=.95  the  pre- 
sent inviscid  calculations  indicate  the  shock  is  at  the  air- 
foil trailing  edge.  It  should  be  noted  that  the  experimental 
data  of  Tijdeman  et  allf*,19  indicate  that  at  Mco=.9  the  shock 
is  between  5 and  10  percent  chord  forward  of  the  aileron 
hinge  point  and  that  the  shock  does  not  reach  the  trailing 
edge  until  very  close  to  Mw=l.  This  tendency  of  the  inviscid 
small  disturbance  result  to  overpredict  the  shock  strength 
and  its  aft  location  is  typical  and  is  most  likely  due  to 
viscous  effects  neglected  in  the  theory.  It  is  recalled  that 
the  comparison  of  the  present  results  with  the  data  (Figure 
36  of  Section  5.3)  is  very  good  for  a fully  supercritical 
Mach  number  of  M^®. 875. 

It  is  noted  that  the  M„,=  .95  pressure  distribution  is 
very  similar  in  shape  to  the  low  supersonic  results  shown  in 
Figure  42.  In  fact,  the  distributions  of  C^-C*  for,  M<x,=  .95, 
1.1,  1.15,  1.2,  compare  to  within  10%  in  eafah  case.  This  re- 
sult is  the  so  called  "Mach  freeze"  phenomena  for  flows  near 
Mach  1 and  in  fact  indicates  that,  although  the  magnitudes  of 
the  pressures  differ  considerably,  the  local  flow  Mach  number 
over  and  in  the  near  field  of  the  airfoil  most  likely  varies 
by  no  more  than  10%  in  the  freestream  Mach  number  range  Ma>= 

.95  to  1.2.  Since  the  nature  of  unsteady  perturbations  to 
the  mean  flow  are  determined  largely  by  the  flow  Mach  number 
field  about  the  airfoil  the  "Mach  freeze"  phenomena  has  im- 
portant consequences  with  respect  to  the  unsteady  perturba- 
tion results.  It  is  expected,  for  example,  that  the  unsteady 
results  for  Moo  between  .95  and  1.2  will  not  vary  a great  deal. 
This  expectation  is  substantiated  somewhat  by  the  unsteady 
aerodynamic  coefficients  presented  in  the  next  section. 

Results  for  unsteady  pressures  due  to  control  surface 
perturbations  of  the  subsonic  mean  flow  solutions  are  pre- 
sented in  Figures  43  through  49.  The  results  in  each  case 
are  compared  to  the  data  of  Tijdeman  et  al  ' Results  for 

the  jump  in  pressure  coefficient  across  the  airfoil  per  unit 
flap  angle  due  to  a quasi-steady  (k=0)  control  surface  de- 
flection are  given  in  Figures  43,44,45,  and  46  for  freestream 
Mach  numbers  of  Mro=.8,  .85,  .9,  .95  respectively.  The  results 
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demonstrate  the  strong  non-uniform  flow  effect  for  transonic 
Mach  numbers.  In  each  case  the  comparison  with  the  data  is 
qualitatively  good  but  quantitatively  in  error.  As  discussed 
above  in  Section  5.2,  the  most  probable  reason  for  the  dis- 
crepancy is  viscous  effects.  Of  some  note  are  the  super- 
critical results  for  M«,=  .85  and  above.  The  results  for  Moo- 
85  (Figure  44)  for  example  shows  that  the  theory  predicts 
the  location  of  the  peak  pressure  perturbation  through  the 
weak  shock  but  overpredicts  its  magnitude.  This  overpredic- 
tion  of  the  perturbation  over  the  non-deflecting  forward 
portion  of  the  airfoil  is  shown  in  this  case  as  well  as  the 
higher  Mach  number  cases  M<»=.9  (Figure  45)  and.Moo-.95  (Fig 
ure  46) . This  indicates  that  in  addition  to  viscous  effects 
wall  interference  effects  could  be  reducing  the  upstream  in- 
fluence in  the  experiment  as  compared  to  the  present  free 
air  results.  The  wind  tunnel  wall  results  of  Section  6.1.4 
add  some  weight  to  this  possibility. 

One  final  comment  with  respect  to  the  quasi-steady 
results  concerns  the  significant  effect  that  the  mean  flow 
results  can  have  on  the  unsteady  perturbations.  A large 
part  of  the  discrepancy  between  the  theory  and  experiment 
in  the  M^.9  and  .95  cases  can  be  attributed  to  the  fact 
that  the  present  steady  theory  predicts  a shock  location 
aft  of  the  experimental  data.  For  example  at  M®-*9  the 
theory  predicts  the  shock  at  the  hinge  point  which  results 
in  a peak  in  the  perturbation  pressure  at  the  hinge  point 
whereas  the  data  shows  the  shock  location  to  be  at  about 
x/c-.l  hence  the  exoerimental  peak  in  perturbation  pressure 
at  that  location.  Also  at  M.-.95  the  present  results  for 
perturbation  pressure  (Figure  46)  compare  very  well,  to  the 
experimental  results  at  M-«1.0.  It  is  recalled  that  the 
theoretical  steady  solution  at  M =.95  compared  more  closely 
(as  far  as  shock  location)  to  the  experimental  result  at 
Moo=l . 0 than  the  Moo=.96  data  (shown  in  the  figure)  for  example. 

Results  for  the  perturbation  pressure  for  fully  un- 
steady, k=. 1 , aileron  oscillations  are  presented  in  Figures 
47,48,  and  49  for  M =.8,  .85,  and  .9  respectively.  The  com- 
parison to  the  data,  shown  in  the  figures,  is  similar  to  the 
comparison  for  the  k=0  results  just  discussed  in  detail.  As 
discussed  there  the  results  are  qualitatively  good  but  quan 
titatively  in  error. 

The  results  presented  in  this  section  indicate  that 
the  most  likely  cause  for  the  discrepancy  between  the  present 
unsteady  perturbation  theory  and  the  data  is  the  effect  ot 
viscosity.  Viscosity  can  have  a direct  effect  on  the  unsteady 
perturbations  but  for  small  perturbations  has  a more  likely 
indirect  effect  by  way  of  its  modification  to  the  mean  flo  . 

In  other  words,  it  is  believed  that  the  present  unsteady 
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results  would  be  in  somewhat  better  agreement  with  the  data 
if  the  steady  results  for  example  for  shock  strength  and  lo- 
cation agreed  with  the  data.  This  brings  up  the  interesting 
possibility  of  improving  the  overall  theory  by  performing 
viscous  modifications  to  the  steady  solution  by  coupling  the 
present  inviscid  theory  to  a boundary  layer  analysis.  Al- 
though this  would  be  a significant  undertaking  it  would  be 
considerably  easier  than  performing  an  unsteady  boundary  layer 
analysis  coupled  to  the  fully  unsteady  theory. 


1 


L 
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Figure  *15.  Airfoil  Pressure  Perturbation  Due  to 
Aileron  Distribution 
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Figure  48.  Airfoil  Pressure  Perturbation  Due  to 


Aileron  Oscillation 
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NACA  64A006 
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Figure  09.  Airfoil  pressure  perturbation  Due  to 
Aileron  Oscillation 
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6.1.2  Summary  of  Unsteady  Aerodynamic  Coefficients 


To  perform  the  transonic  flutter  analysis  of  Section  7.0, 
unsteady  aerodynamic  coefficients  for  pitch,  plunge  and  control 
surface  oscillations  have  had  to  be  generated  for  the  Mach 
number  range  Mn  3 .8  to  1.2  and  low  reduced  frequency  range 
(k  = 0.0  to  0.2).  The  results  were  calculated  as  unsteady  per- 
turbations to  the  zero  angle  of  attack  steady  solutions  pre- 
sented above.  The  resulting  unsteady  aerodynamic  coefficients 
for  each  degree  of  freedom  are  presented  in  this  section  in 
graphical  form  as  functions  of  free stream  Mach  number  for  each 
reduced  frequency  considered.  The  definitions  used  for  the 
coefficients  are  believed  consistent  with  the  American  litera- 
ture and  are  defined  in  detail  in  Appendix  B,  which  also  con- 
tains tables  of  the  numerical  results. 

The  entire  matrix  of  cases  for  this  airfoil  involved 
approximately  110  separate  and  distinct  steady  and  unsteady 
calculations.  This  matrix  of  calculations  required,  by 
conservative  estimate,  approximately  1.25  hours  of  CDC  7600 
time.  This  is  a better  measure  of  the  efficiency  and  utility 
of  the  method  than  quoting  the  often  widely  varying  run  time 
for  an  individual  case.  It  also  indicates  that  transonic 
flutter  analyses  for  two-dimensional  airfoils  are  possible  with 
relatively  modest  expenditures  of  computer  time. 

Figures  50  through  59  present  a summary  of  the  results 
for  real  and  imaginary  parts  of  the  lift,  moment  and  hinge 
moment  coefficients  per  unit  flap  angle,  pitch  angle,  or  plunge 
deflection  for  the  three  rigid  body  degrees  of  freedom.  The 
figures  are  presented  for  completeness  and  the  the  readers 
perusal  in  total  as  they  point  out  the  important  transonic 
effects  on  all  coefficients.  Figure  50,  for  example,  presents 
the  lift  coefficient  due  to  control  surface  oscillation  and 
demonstrates  the  expected  strong  nonlinear  effects  in  the  trans- 
onic regime  as  well  as  the  large  phase  shifts  in  unsteady  forces 
Portions  of  the  curves  between  Mo,  - .95  and  1.1  are  extrapolated 
due  to  the  inability  to  perform  stable  unsteady  calculations  in 
this  region.  For  example  at  M«,  = .95  the  numerical  method  was 
unstable  for  k > .1  and  for  all  reduced  frequencies  at  Moo  = 1.05 
The  required  extrapolations  are  made  with  some  degree  of  con- 
fidence however  because  of  the  "Mach  freeze"  effect  in  this 
region.  As  discussed  just  above,  the  steady  calculations  for 
Moo  = .95  and  1.05  show  only  slight  differences  in  the  flow  Mach 
number  in  the  solution  field  on  or  near  the  airfoil  so  that 
unsteady  perturbations  are  expected  to  vary  only  slightly  from 
the  Mo,  = .95  or  1.1  results.  This  should  be  expecially  true 
for  the  control  surface  perturbations  as  in  fact  is  verified 
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by  Figures  50,  51  and  52.  Figures  53  through  59  show  that  the 
variation  in  the  coefficients  between  Mffi  = .95  and  1.01  for  the 
pitch  and  plunge  nodes  is,  as  night  be  expected,  sonewhat  greater 


ICIENT  PER  UNIT  FL£PS  ANGLE 
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6.1.3  Mean  Angle  of  Attack  Effects  on  Unsteady  Aerodynamic 
Coefficients 


The  basic  theory  and  all  of  the  unsteady  results  pre- 
sented  in  this  report  indicate  and  demonstrate  that  non-uniform 
mean  flow  effects  dominate  the  development  of  the  unsteady 
flow  at  transonic  speeds  even  in  the  limit  of  small  unsteady 
perturbations.  This  is  in  contradistinction  to  the  usual 
experience  for  subsonic  and  supersonic  flows  for  which  the 
first  order  unsteady  flow  uncouples  and  is  therefore  inde- 
pendent of  the  mean  flow.  One  manifestation  of  the  non-uniform 
f low  effects  is  the  effect  of  mean  airfoil  angle  of  attack 
(or  mean  flap  angle)  which  since  it  may  significantly  affect 
the  mean  flow  can  therefore  influence  the  resulting  unsteady 
perturbation.  This  effect  is  briefly  demonstrated  here  for 
the  NACA  64A006  at  a supercritical  Mach  number  of  .ob. 

Figure  60  below  presents  the  calculated  pressure  distri- 
bution for  this  airfoil  at  = .35  for  mean  angles  of  attack 
a = o,  1°  and  2°.  The  figure  shows  that  the  upper  and  lower 
surface  pressure  distributions  which  are  identical  for  aQ  - 0 
(symmetric  airfoil)  are  significantly  altered  for  angles  of 
attack  of  1°  and  2°.  As  shown,  as  angle  of  attack  increases, 
the  lower  surface  becomes  subcritical  and  the  upper  surface 
becomes  more  and  more  supercritical.  Also  the  relative  y 
weak  shock  at  an  = 0 becomes  stronger  and  moves  aft  on  tne 
upper  surface  much  as  it  would  if  the  freestream  Mach  number 
were  increased.  It  seems  clear  that  this  relatively  large 
change  in  the  mean  flow  will  significantly  affect  the  unsteady 

flow  perturbation. 

This  is  indeed  shown  to  be  the  case  in  Figure  61,  62, 
and  63  which  present  the  variation  with  reduced  frequency  of 
the  lift,  moment  and  hinge  moment  coefficients,  respectively, 
due  to  the  oscillating  aileron.  Results  are  presented  there 
for  0,  1°  and  2°  mean  angle  of  attack  and  limited  experimental 
data  at  zero  angle  of  attack  are  presented  for  comparison.  The 
strong  reduced  frequency  effect  shown  for  and  ? g 
weak  dependence  for  Cng  are  consistent  with  Subsonid  linear 
theory.  The  primary  result  of  note  is  the  prediction  of  a 
substantial  mean  angle  of  attack  effect.  Of  some  considerable 
interest  to  the  flutter  problem  is  the  change  in  sign  of  the 
imaginary  (out  of  phase)  part  of  the  hinge  moment  due  to 
mean  angle  of  attack.  This  indicates  a potential  for  one 
degree  of  freedom  flutter  (control  surface  buzz  due  solely  to 
mean  airfoil  angle  of  attack;  an  important  result  not  predicted 
by  linear  theory. 
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NACA  64A006 


Figure  60,  Steady  Pressure  Coefficients 
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Figure  61.  Mean  Angle  of  Attack  and  Reduced  Fre 
quency  Effect  on  Lift  Coefficient 
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Figure  62,  Mean  Angle  of  Attack  and  Reduced  Fre- 
quency Effect  on  Moment  Coefficient 
(About  Leading  Edge) 
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6.1.4  Wall  Interference  Effects  on  Unsteady  Aerodynamic 
Coefficients 


The  extentions  of  the  steady  and  unsteady  flow  solu- 
tion method  to  include  wall  effects  were  described  above  in 
Section  3.3.  In  this  section  the  capability  is  used  to  brief- 
ly examine  the  potential  effect  of  the  wind  tunnel  configura- 
tion used  in  the  Tijdeman  et  al.  experiments  on  their  results. 

The  approximate  modeling  of  the  wall  boundary  condition 
was  discussed  in  the  earlier  section  and  steady  and  quasi- 
steady perturbation  results  for  the  NACA  64A006  were  presented 
there.  These  results  indicated  that  although  the  effect  of  the 
walls  was  appreciable  there  was  little  difference  between  solid 
and  the  slotted  wall  results.  Also  as  shown  there  the  inclu- 
sion of  wall  effects  seemed  in  fact  to  increase  the  discrepancy 
with  experiment  with  respect  to  the  important  peak  pressure 
perturbation  through  the  shock  (re:  Figure  18).  This  is  shown 

in  Figure  64  below  which  essentially  repeats  Figure  18  with  the 
addition  of  experimental  data  for  the  same  airfoil  and  a quasi- 
steady aileron  deflection  at  Mo  = .875.  As  shown  in  Figure  6 
the  addition  of  the  slotted  walls  reduces  the  agreement  with 
the  data  with  respect  to  location  and  magnitude  of  the  peak 
pressure  through  the  shock  but  enhances  the  agreement  with  res- 
pect  to  the  depressed  pressure  perturbation  forward  of  the  peaK. 
The  reasons  for  the  difference  between  the  free  air  and  the  wall 
solution  can  be  explained  almost  entirely  by  the  effect  of  the 
walls  on  the  corresponding  steady  solution.  That  is,  the  effect 
of  the  wall  is  to  strengthen  the  shock  and  move  it  aft  on  the 
airfoil  hence  the  enhanced  perturbation  pressure  through  the 
shock  and  the  aft  location  with  respect  to  the  free  air  result. 
The  depressed  perturbation  forward  of  the  peak  is  due  to  the 
stronger  shock  and  larger  supersonic  region  in  the  steady  wall 
solution,  both  of  which  have  the  effect  of  decreasing  the  sig- 
nal which  reaches  the  forward  part  of  the  airfoil  from  the  con- 
trol surface. 


The  data  at  Ma  = .875  is  added  to  present  a possible  ex- 
planation for  the  discrepancy  between  the  present  results  and 
the  data.  Although  it  may  be  fortuitous,  the  comparison  of  the 
present  porous  wall  results  at  Moo  = .85  and  the  data  at  M™  .875 
very  good  with  the  exception  of  an  overprediction  of  the  magni- 
tude  of  the  perturbation  through  and  just  aft  of  the  shock.  As 
discussed  earlier,  the  inviscid  small  disturbance  steady  results 
are  almost  without  exception  more  supercritical  than  the  corres- 
ponding experimental  results.  This  is  usually  attributed  to  the 
effect  of  viscosity,  neglected  in  the  theory,  which  has  the 
effect  of  weakening  and  moving  the  inviscid  shock  forward  on  tne 


is 
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airfoil.  If  this  one  comparison  is  to  be  believed,  a possible 
alternative  to  the  viscous  coupling  suggested  above  would  be 
to  include  wall  effects  in  the  inviscid  theory  and  to  account 
for  viscosity  as  an  effective  "shift”  in  the  Mach  number  for 
comparison  to  experimental  data.  This  could  clearly  only  be 
done  apriori  and  of  course  would  reduce  the  practical  utility 
of  the  method. 

Figures  65,  66,  and  67  are  presented  to  show  the  effect 
of  solid  or  ventilated  walls  on  the  lift  coefficient  due  to 
control  surface,  pitch  and  plunge  oscillations  respectively. 

As  shown,  the  effect  is  quite  appreciable  in  the  low  reduced 
frequency  regime.  Of  some  particular  note  is  the  fact  that 
the  trend  with  respect  to  reduced  frequency  is  different  for 
the  wall  cases  as  compared  to  the  free  air  case.  This  could 
be  due  to  the  fact  that  the  wall  cases  are  effectively  for  a 
higher  Mach  number  due  to  the  influence  of  the  walls  on  the 
steady  flow. 


Figure  64.  Effect  of  Wind  Tunnel  Walls  on  Airfoil 
Perturbation  Pressure 
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Figure  65.  Wind  Tunnel  Wall  and  Reduced  Fre- 
quency Effect  on  Lift  Coefficient 
Due  to  Control  Surface  Rotation 
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6.2 


Results  for  a NACA64A410  Airfoil 


Steady  and  unsteady  results  are  presented  in  this  sec- 
tion for  a NACA64A410  airfoil  for  supercritical  Mach  numbers. 
Calculations  have  been  performed  for  pitch  and  plunge  oscill- 
ations at  reduced  frequencies  of  k=0  to  0.2  for  free  air  Mach 
numbers  of  M^.75,  .8  and  .9  at  zero  mean  angle  of  attack  and 
additionally  at  1°  and  2°  angles  of  attack  at  Moo=.8.  The  sol- 
utions for  this  airfoil  follow  the  same  general  trends  as 
those  already  discussed  in  detail  for  the  NACA64A006  so  that 
the  results  are  presented  for  completeness  with  minimal  addi- 
tional comment. 

The  NACA64A410  airfoil  has  a 10%  thickness  and  adds 
the  additional  variable  of  a 4%  camber,  both  of  which  signif- 
icantly increase  the  non-uniform  flow  effects  in  comparison 
to  the  64A006  which  is  of  the  same  family  of  airfoils.  The 
steady  results  are  shown  in  Figures  68  and  69  below  and  dem- 
onstrate the  large  supercritical  region  even  at  a relatively 
low  Mach  number  of  M^.75.  Due  to  the  airfoils  non-symetri- 
cal  shape,  there  is  a relatively  large  difference  between  the 
upper  and  lower  surfaces  at  zero  angle  of  attack  as  shown  in 
Figure  68.  The  results  for  progressively  higher  Mach  numbers 
show  the  strengthening  and  aft  movement  of  the  shock  and  it 
is  noted  that  the  shock  is  already  at  the  airfoil  trailing 
edge  at  M<x>=.85.  The  results  for  different  angles  of  attack, 
given  in  Figure  69,  indicate  that  at  least  for  small  angles 
of  attack  the  variation  of  the  steady  solution  is  also  small 
due  presumably  to  the  fact  that  the  airfoil  is  already  at  a 
relatively  large  "effective"  angle  of  attack  because  of  its 
camber.  Based  on  these  comparisons  it  is  expected  that  the 
unsteady  perturbations  to  these  mean  flows  would  show  a 
large  Mach  number  variation  but  perhaps  a small  effect  of 
angle  of  attack.  A table  of  the  calculated  steady  lift  co- 
efficients is  presented  below  for  all  the  cases  calculated. 

The  values  demonstrate  the  large  transonic  flow  effect  on 
lift.  In  particular  the  zero  angle  of  attack  lift  is  shown  to 
peak  at  about  M =.85  and  to  decrease  significantly  at  M^.9. 
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.918 

1.047 

.773 

.305 

Table  2.  Lift  Coefficients  for  NACA  64A410 


•phe  results  for  unsteady  lift  and  moment  coefficients 
due  to  pitch  and  plunge  perturbations  to  the  above  steady 
solutions  are  presented  in  Figures  70  through  73.  The  re- 
duced frequency  variation  of  the  various  coefficients  as 
shown  in  the  figures,  is  typical  of  those  shown  previously 
and  again  demonstrates  the  large  reduced  frequency  effect 
in  the  low  frequency  regime.  Also  of  considerable  note  is 
the  large  non-linear  effect  of  mean  flow  Mach  number.  Another 
result  not  shown  is  that  the  efffect  of  mean  angle  of  attack 
at  M =.8  is  quite  small.  In  fact  the  aerodynamic  coefficients 
vary°°by  only  a few  precent  for  angles  of  attack  between  0 
and  2°. 
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Mach  Number  and  Reduced  Frequency 
Effect  on  Moment  Coefficient 
(About  X=0)  Due  to  Pitch 


Figure  71 
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Figure  72.  Mach  Number  and  Reduced  Frequency 
Effect  on  Lift  Coefficient  Due  to 
Plunge 


7.0  TRANSONIC  AIRFOIL  FLUTTER  ANALYSIS 


The  primary  purpose  of  the  research  reported  here  was 
to  develop  a method  for  calculating  unsteady  aerodynamic  coef 
ficients  for  use  in  transonic  flutter  analyses.  These  analy- 
ses were  carried  out  for  two  representative  airfoils  to  illus 
trate  the  practical  effects  of  transonic  unsteady  aerodynamic 
on  flutter  computations.  The  flutter  analysis  method  and  re- 
sults are  summarized  in  this  section 


mean 

position 


mid  chord 

elastic  axis 


center  of  mass  ** 

hinge  point 


Figure  74,  Notation  for  Flutter  Analysis 


37 

The  flutter  analysis  follows  the  method  of  Fung  and 
is  the  so-called  V-g  method  which  is  commonly  used  in  current 
practice.  A two-dimensional  airfoil  with  semi-chord,  b,  and 
of  unit  length  in  the  spanwise  direction  is  considered  free 
to  oscillate  in  three  degrees  of  freedom  (plunge,  pitch,  and 
control  surface  oscillation)  as  shown  schematically  in  Fig- 
ure 74  where: 


h = plunge  deflection  of  the  elastic  axis,  positive 
downward 

a = pitch  about  the  elastic  axis,  positive  nose  up 


angular  deflection  of  aileron  about  the  aileron 
hinge  point  relative  to  the  wing  chord,  positive 
tail  down. 


The  notation  used  in  the  analysis,  some  of  which  is  defined 
in  the  figure,  is  consistent  with  that  used  by  Fung.  It  is 
noted  that  all  dimensions  are  non-dimensionalized  with  respect 
to  the  semi-chord  and  distances  are  measured  with  the  mid-chord 
point  as  the  origin. 

It  is  assumed  that  each  degree  of  freedom  is  restrained 
by  a linear  elastic  restraint  with  zero  structural  damping. 

The  three  equations  of  motion  are  determined  by  summing  the 
inertia,  elastic  and  aerodynamic  forces;  and  assuming  harmonic 
oscillations  these  equations  at  the  flutter  point  can  be  writ- 
ten in  the  following  matrix  form: 


/ \ ( 

h)  (h) 

l\i  k2  [M]  - [A]  J ^ 

a > = UK]  < 

8 ) ( B) 

(52) 


where  the  flutter  eigenvalue  is: 


X = m (l+ig)ui2  b2/U: 


(53) 


with  flow  velocity  and  density,  U and  p respectively,  refer- 
ence frequency,  uir,  wing  mass  per  unit  span,  m and  mass  ra- 
tio  v = m/71  p b!.  Also,  the  artificial  structural  damping 
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required  to  maintain  simple  harmonic  motion  is  denoted  by  g 
in  the  flutter  eigenvalue.  The  reduced  frequency  of  oscilla- 
tion is  defined  in  the  section  using  the  usual  American  defi- 
nition, namely: 


k = 


tub 

U 


(54) 


where  it  is  noted  that  the  semi-chord,  b,  is  used  in  the  ex- 
pression rather  than  the  chord,  c,  used  in  other  sections  of 
this  report.  It  is  hoped  that  this  change  of  definitions  for 
k does  not  cause  undue  confusion.  The  mass  matrix  [M]  and 
stiffness  matrix  [K]  are  written  in  terms  of  the  airfoil  in- 
ertial and  stiffness  parameters  as  follows: 


(55) 


and 


(K] 


(57) 


where  many  parameters  have  been  defined  in  Figure  74  and  where 
Xo  is  the  distance  from  the  aileron  hinge  point  to  its  center 
or  mass  and  ra  and  rg  are  the  radii  of  gyration  for  the  mass 
moment  of  inertia  of  the  airfoil  about  the  elastic  axis  and 
for  the  aileron  about  its  hinge.  The  frequencies,  u)a,  u>g 
are  the  uncoupled  natural  frequencies  for  the  airfoil  in  plunge, 
pitch  and  control  surface  oscillation  respectively.  Finally, 
the  aerodynamic  coefficient  matrix  is  given  by: 


which  transforms  the  unsteady  aerodynamic  coefficients  as  de- 
fined in  this  report  (see  Appendix  B)  to  account  for  the  non- 
dimensionalization  with  respect  to  b used  in  this  section 
(rather  than  c used  elsewhere  in  the  report) . The  transfor- 
mation also  accounts  for  the  shift  in  the  definition  of  air- 
foil moments  from  the  leading  edge  (as  used  in  the  definition 
of  Cm  , etc.)  to  the  elastic  axis,  ah#  required  in  this  analy- 
sis. a 

Performing  a standard  eigenvalue  analysis  of  the  above 
system  (Eq.  52) , provides  the  formalism  for  determining  the 
so-called  flutter  point.  For  a given  airfoil,  the  flutter 
point  at  fixed  Mach  number  and  air  density  is  determined  by 
finding  the  lowest  value  of  flow  velocity  for  which  the  arti- 
ficial damping  parameter,  g (imaginary  part  of  X)  is  zero. 

This  is  found  by  solving  for  the  eigenvalues  of  Eq.  52  for 
various  reduced  frequencies  and  cutting  and  trying  until  the 
flutter  velocity  is  determined  within  acceptable  accuracy. 

The  results  of  this  analysis  are  presented  here  in  the 
form  of  a dimensionless  flutter  speed,  u*,  and  flutter  fre- 
quency, <jj*,  where 


u* 


U 


b a). 


a 


(j 


* = 


u> 


a 


(58) 


Representative  values  for  the  airfoil  geometry,  inertia  and 
stiffness  parameters  were  used  as  follows: 


-134- 


Variations  from  these  values  were  used  to  illustrate  special 
trends  as  will  be  discussed. 


The  primary  results  of  the  analysis  are  for  the  two 
degree  of  freedom  (2  d.o.f.)  case  for  the  NACA  64A006  and 
these  results  are  presented  in  Figure  75  through  77.  In  the 
2 d.o.f.  case  the  aileron  is  assumed  fixed  with  respect  to 
the  airfoil  so  that  the  motion  is  restricted  to  pitch  and 
plunge  oscillations.  Figures  75  and  76  present  results  for 
flutter  velocity  and  frequency  respectively,  as  functions  of 
freestream  Mach  number.  The  figures  demonstrate  the  crucial 
importance  of  the  transonic  Mach  number  range  by  way  of  the 
dip  in  flutter  velocity  at  M^  = .85  and  the  corresponding 
large  variation  in  flutter  frequency  in  this  range.  The  re- 
sults for  flutter  velocity  are  presented  for  two  values  of  the 
center  of  mass  location,  x , which  indicate  the  importance  of 
this  parameter  on  flutter  speed.  It  is  noted  that  the  airfoil 
is  statically  stable  (divergence  free)  for  Mach  numbers  great- 
er than  Mo,  - .85  due  to  the  shift  in  the  aerodynamic  center 
aft  of  the  elastic  axis  for  M*  > .85.  Also  note  that  the 
normal  flutter  mode  for  all  cases  given  in  the  figures  was 
dominated  by  motion  in  plunge  for  wh/wa  = .2.  It  would  be 
expected  that  the  coupling  between  the  pitch  and  plunge  mode 
would  increase  as  uih  approaches  w0  and  that  therefore  the  cri- 
tical flutter  velocity  would  be  decreased.  This  effect  is 
shown  in  Figure  77  for  wh/u)  up  to  0.5. 


The  effect  of  mean  angle  of  attack  and  wind  tunnel 
wall  interference  on  the  critical  flutter  parameters  was 
also  examined  for  the  limited  cases  for  which  the  aerodynam- 
ic coefficients  have  been  calculated.  The  results  are  sum- 
marized in  Table  3 for  the  NACA64A006  at  Mo<>=.85  and  in- 
dicate a modest  beneficial  effect  of  mean  angle  of  attack  on 
flutter  velocity  and  a somewhat  greater  ■'detrimental  effect 
of  wall  interference.  Both,  however,  have  a considerably 
greater  effect  on  the  flutter  frequency  which  is  decreased 
from  w*  = .473  for  the  free  air  case  to  w*  = .075  for  the 
wall  case. 

Fully  three  degree  of  freedom  flutter  calculations 
have  also  been  performed  for  the  NACA  64A006  with  not  alto- 
gether complete  results.  Due  to  the  rather  restricted  re- 
duced frequency  range  (0  < k < .1)  of  the  present  aerodynam- 
ic data,  it  was  impossible  to  determine  the  true  flutter 
velocity  because  an  aileron  dominated  normal  mode  (aileron 
buzz)  was  fluttering  at  a frequency  which  was  out  of  the 
range  of  the  data.  It  was  determined,  however,  that  the 
aileron  significantly  damped  the  plunge  dominated  normal 
mode  which  fluttered  in  the  2 d.o.f.  case.  This  happened 
even  for  unrealistically  high  values  of  u)g/wa  ^ 50  and  for 
values  of  wg/wa  30  the  plunge  dominated  mode  was  complete- 
ly damped  and  therefore  did  not  flutter.  This  curious  coup- 
ling of  the  aileron  and  plunge  modes  and  the  high  frequency 
aileron  buzz  phenomena  would  seem  to  be  peculiar  to  the  pre- 
sent unsteady  aerodynamic  data  and  merits  further  investiga- 
tion. Time  constraints  on  the  present  study  did  not  permit 
a comprehensive  study  of  the  3 d.o.f.  flutter  case. 

Finally,  it  is  noted  that  2 d.o.f.  flutter  calculations 
have  been  performed  for  the  NACA  64A410.  However,  the  re- 
stricted Mach  number  range  of  the  aerodynamic  data  did  not  per- 
mit an  in-depth  examination  of  the  transonic  flow  effects. 

For  example,  at  M«  = .75  the  airfoil  flutter  velocity  was  25.1 
and  flutter  frequency  of  .19.  Mean  angle  of  attack  had  the 
expected  negligible  effect  on  the  flutter  parameters  as  did 
the  data  at  Mo,  = .8  and  .85. 

In  summary,  the  two  degree  of  freedom  flutter  results 
for  the  NACA  64A006  airfoil,  presented  here,  demonstrate  the 
crucial  importance  of  the  transonic  speed  range  with  respect 
to  airfoil  flutter.  A thorough  examination  of  three  degree 
of  freedom  flutter  would  require  the  generation  of  unsteady 
aerodynamic  coefficients  over  a more  extended  range  of  reduced 
frequencies  than  were  considered  in  the  present  work. 
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Figure  75.  Flutter  Velocity  Versus  Mach  Number 
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Figure  76.  Flutter  Frequencies  versus  Mach  Number 


8 , 0 CONCLUSIONS 


' In  this  report,  a treatment  of  unsteady  transonic  flow 

ias  a small  perturbation  about  the  nonlinear  small  disturbance 
steady  flow  is  presented.  The  technique  has  been  applied  to 
representative  airfoils  oscillating  at  subsonic,  through  tran 
sonic  to  supersonic  Mach  numbers.  It  is  shown  that  the  present 
theory  matches  linearized  unsteady  subsonic  or  supersonic  theory 
as  Mach  number  decreases  or  increases  away  from  sonic  respec- 
tively. The  theory  includes  the  effects  of  thickness  and  angle 
of  attack  which  are  ignored  in  linear  theories  but  which,  as 
the  results  presented  here  demonstrate,  are  significant  in  the 
low  frequency  unsteady  transonic  speed  regime.  A measure  of 
the  efficiency  of  the  computational  scheme  is  given  by  the  fact 
that  the  matrix  of  calculations  required  to  generate  unsteady 
aerodynamic  coefficients  for  a transonic  flutter  study  of  a 
three  degree  of  freedom,  two  dimensional  airfoil  can  be  com- 
pleted in  about  1 and  1/4  hours  of  CDC  7600  time. 

As  a result  of  the  research  reported  here,  the  method 
has  been  generalized  to  treat  higher  frequency  effects,  super- 
sonic freestream  flows,  wind  tunnel  walls  and  to  perform  three 
dimensional  planar  wing  calculations.  The  method  is  currently 
in  a state  of  development  that  the  study  of  such  effects  is 
routine  using  computer  programs  developed  during  the  study. 

The  method  is  somewhat  limited  by  a numerical  instability 
which  develops  for  freestream  Mach  numbers  very  near  Mach  1, 
but  a wide  range  of  supercritical  flows  may  be  successfully 
calculated.  A possible  resolution  to  this  instability  has 
been  postulated  based  on  an  inner/outer  expansion  analysis 
for  near  sonic  Mach  numbers  but  its  implementation  would  re- 
quire further  development  effort. 

The  results  presented  here  indicate  that  the  theory 
and  numerical  solution  method  provide  a meaningful  represent- 
ation of  inviscid  transonic  flows  about  practical  airfoils 
for  reasonable  amplitudes  of  unsteady  motion.  The  results 
are  shown  to  compare  adequately  with  the  "exact"  numerical 
calculations  or  Magnus  and  Yoshihara  for  Low  frequency  un- 
steady supercritical  flows,  for  an  expenditure  of  over  two 
orders  of  magnitude  less  computer  time.  Although  these  com- 
parisons show  some  qualitative  differences  in  unsteady  airfoil 
pressure  perturbations,  unsteady  aerodynamic  forces  compare 
very  well  in  both  amplitude  and  phase.  Detailed  comparisons 
of  the  present  results  to  limited  available  two  dimensional 
experimental  data  shows  good  qualitative  agreement.  Quanti- 
tative discrepancies  would  seem  to  be  related  mainly  to  viscous 
effects,  not  accounted  for  in  the  model,  and  to  possible  wind 
tunnel  wall  or  three  dimensional  flow  effects  in  the  data. 
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A brief  study  of  such  effects  indicates  that,  as  expected, 
they  can  be  quite  significant  in  the  transonic  speed  range. 

The  need  for  an  accurate  and  efficient  method  for  cal- 
culating unsteady  aerodynamic  coefficients  at  transonic  speeds 
is  pointed  out  by  the  flutter  dynamics  study  presented  in  this 
report.  Aerodynamic  coefficients  were  calculated  using 
present  method  and  used  in  a flutter  analysis  of  a JJACA  64A 
with  three  rigid  body  degrees  of  freedom  (pitch,  plunge, 
trll  surface dotation)  and  of  a NACA  64A410  with  two  rigid  body 
degrees  of  freedom  (pitch,  plunge)  The  results: speed 
pected  critical  effect  of  the  transonic  regime  on  flutter  speed. 

It  seems  clear  that  viscous  effects  are  quite  important 
for  unsteady  flows  in  the  transonic  speed  regime.  The  compar- 
ison of  the  present  inviscid  results  to  data,  given  above,  in 
dieate  that  the  inclusion  of  viscous  effects  is  crucial  to  the 
practical  accuracy  of  an  unsteady  transonic  theory.  The  modest 
computer  time  requirements  of  the  present  small  disturbance 
theorymake  it  a^otential  candidate  for  the  introduction  of 
such  effects  and  it  is  recommended  that  this  be  the  next 
for  further  development  of  the  method. 
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APPENDIX  A 


FARFIELD  BOUNDARY  PRESCRIPTIONS 


In  the  solution  technique  presented  in  this  report, 
a solution  field  of  finite  extent  is  defined.  The  comple- 
tion of  the  finite  difference  analogs  to  the  steady  and  un- 
steady transonic  flow  problems,  require  that  boundary  con- 
ditions be  fixed  in  some  manner  on  the  farfield  boundary. 

In  this  report,  calculations  are  performed  for  two  dimen- 
sional airfoil  sections  in  subsonic  or  supersonic  free  air 
or  constrained  by  solid  or  porous  wind  tunnel  walls  and  for 
three  dimensional  planar  wings  in  subsonic  free  air  flows. 
In  each  case  a different  farfield  prescription  is  required 
and  the  purpose  of  this  appendix  is  to  summarize,  for  com- 
pleteness sake,  the  asymptotic  solutions  used  to  fix  these 
farfield  boundary  conditions. 


A.l  Two-Dimensional  Airfoil  in  Free  Air 


The  basic  method  used  to  develop  asymptotic  solutions 
for  the  steady  or  unsteady  perturbation  potentials  involves 
the  integral  representation  which  results  from  applying 
Green's  theorem  to  the  appropriate  partial  differential 
equation.  Suitable  approximations  are  made  to  the  integral 
expressions  valid  far  from  the  airfoil  and  a farfield  solu- 
tion results.  Results  for  steady  and  unsteady  free  air  flows 
are  presented  below. 


A. 1.1  Steady  Free  Air  Flows 


The  farfield  for  the  non-linear  steady  potential  was 
examined  in  the  manner  mentioned  above  initially  by  Cole-32 
and  most  extensively  by  Klunker29f  whose  result  is  used. 

The  farfield  representation  describes  thickness  and  lift  ef- 
fects to  leading  order  and  is  given  by: 
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where  t(£)  is  the  airfoil  thickness  distribution.  The  doublet 
strength  due  to  the  nonuniform  flow  is  given  by  the  integral  of 
(<f>x) 2 over  the  solution  field.  This  is  evaluated  periodically 
as  the  iterative  solution  pxoceeds  using  a trapezoidal  rule  in- 
tegration scheme.  Similarly# in  calculations  with  lift,  the  air- 
foil circulation#  y,  is  updated.  As  the  doublet  strength  and 
circulation  are  refined#  Eq.  A.l  is  used  to  update  <J>°  on  the 
boundaries  of  the  finite  difference  grid.  Thus  the  farfield 
boundary  condition  depends  on  the  numerical  solution  and  is 
treated  as  a Dirichlet  boundary  condition. 

For  supersonic  freestream  Mach  numbers#  an  approximate 
characteristic  relation  can  be  developed  from  either  the  inte- 
gral equation  or  the  governing  partial  differential  equation. 
This  was  done  for  steady  flows  by  Murman22  and  for  unsteady 
flows  by  the  present  authors (2°) . This  characteristic  condi- 
tion is  thus  used  to  define  a boundary  condition  on  the  farfield 
boundaries.  It  is  assumed  that  the  disturbance  from  the  body  is 
weak  so  that  the  characteristic  relation  on  the  incoming  charac- 
teristic applies.  The  first  order  approximation  to  the  relation 
is  given  and  the  farfield  boundaries  are  treated  as  follows: 


<t>°.  = + /-K  <|>°  on  y = ± °° 
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<j>°  = = 0 as  x 


(A. 2) 
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as  x -*■  +°° 


The  upstream  boundary  (x  •>  -®>)  is  taken  ahead  of  the  bow  shock 
so  that  the  flow  is  uniform.  Because  of  the  backward  differen- 
cing used  for  supersonic  grid  points,  no  prescription  is  re- 
quired on  the  downstream  boundary  (x  -*•  +°°)  as  long  as  the 
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boundary  is  completely  supersonic.  If  during  the  iteration 
procedure  a grid  point  on  the  boundary  becomes  subsonic,  a 
zero  gradient  (4>x  - 0)  condition  is  used.  The  above  condi- 
tions are  incorporated  in  the  finite  difference  procedure  as 
a Dirichlet  condition  on  the  upstream  boundary  and  using  one 
sided  finite  difference  equations  for  the  gradient  conditions 
on  the  top  and  downstream  boundaries. 


A. 1.2  Unsteady  Free  Air  Flows 


In  the  same  manner  as  the  steady  flow,  the  farfield 
boundary  conditions  for  the  unsteady  perturbation  potential 
in  the  low  frequency  regime  are  defined  from  an  asymptotic 
expression  for  subsonic  freestream  flow  and  an  approximation 
to  the  unsteady  characteristic  relation  for  supersonic  free- 
stream flows.  For  subsonic  flow  the  asymptotic  expression  is 
used  to  fix  values  of  <f>  in  the  Dirichlet  sense  on  all  four 
boundaries  of  the  grid.  The  expression  was  derived  in  Refer- 
ence 1 and  is  given  by: 
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where  k = J2/K,  R,=  k \/(x-l)z+Ky2  and  Hi  is  the  Hankel  function 
of  the  second  kind  of  ith  order.  The  integral  cannot  be  eval- 
uated in  closed  form  but,  since  it  is  over  a finite  range, 
it  can  be  integrated  numerically  with  little  difficulty.  It 
is  noted  that  the  farfield  approximation  depends  on  the  solu- 
tion through  the  airfoil  circulation  y and  also  through  the 
"airfoil  integral".  Thus  as  the  airfoil  circulation  and  po- 
tential distributions  are  refined,  the  farfield  is  periodi- 
cally updated  using  equation  A.3-A.5. 

Supersonic  freestream  flows  are  treated  in  a analagous 
manner  to  the  steady  flow.  A low  frequency  unsteady  character- 
istic condition  has  been  developed^®  and  is  used  along  with  up- 
stream and  downstream  conditions  as  farfield  boundary  conditions. 

The  relations  are: 
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Two-Dimensional  Airfoil  Between  Wind  Tunnel  Walls 


A.  2 


The  effect  of  wind  tunnel  walls  in  the  steady  and  un- 
steady flow  is  of  considerable  concern  especially  for  tran- 
sonic flows  and  the  effect  can  be  modeled  by  treating  the 
appropriate  boundary  conditions  at  the  wall  and  far  upstream 
and  downstream. 


A linearized  boundary  condition  which  relates  the  nor- 
mal velocity  at  the  wall  to  the  local  pressure  is  used  in  the 
present  study  for  general  porous  walls  (solid  walls  are  con- 
sidered as  a special  case) . Also  asymptotic  expressions  valid 
for  upstream  and  downstream  have  been  developed  for  subsonic 
upstream  flow.  In  the  present  work  a farfield  expression  for 
solid  walls  has  been  developed  and  is  used  in  the  general  case. 
The  expression  results  from  the  superposition  of  a infinite 
sequence  of  doublets  and  vortices  (of  alternating  signs)  each 
separated  by  the  scaled  distance  between  the  walls,  h.  This 
solution  is  used  to  fix  a Dirichlet  boundary  condition  on  the 
upstream  and  downstream  boundaries  of  the  finite  difference 
grid.  The  resulting  set  of  boundary  conditions  for  subsonic 
steady  flow  are: 
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and  where  P is  the  wall  porosity  coefficient.  As  before,  the 
doublet  strength,  V,  and  the  airfoil  circulation  y are  refined 
as  the  solution  proceeds  and  used  to  update  the  farfield.  As 
stated  above,  the  wall  boundary  condition  is  treated  as  the 
first  order  (linear)  relation  between  the  perturbation  velocity 
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normal  to  the  wall  and  the  local  pressure  coefficient  (C  = 
-2<f>x),  assuming  the  plenum  pressure  is  equal  to  the  freeitream 
static  pressure.  For  solid  wall  note  that  P=0  so  that  <f>°=0 
at  the  wall.  An  Approximation  for  P for  general  porous  £r 
slotted  walls  can  be  developed  by  assuming  the  boundary  condi- 

^Sothe  wei9hted  average  of  a free  jet  ($®=0)  and  a solid 
wall  (4>y=0)  where  the  weighting  factors  are  the  open  area 

ratio,  S,  or  1-S  respectively.  The  resulting  expression  for 
* is : 


P 


S 1 

1-S  [(1+y)6m£] 1/3 


(A. 8) 


The  boundary  conditions  for  the  wall  case  with  super- 
sonic upstream  flow  can  be  treated  in  a simple  manner  as  long 
as  the  upstream  and  downstream  boundaries  are  defined  far 
enough  from  the  airfoil.  In  this  case  the  wall  condition  is 
the  same  as  given  in  Equation  A. 8 and  the  upstream  condition 
is  taken  as  uniform  flow.  No  downstream  condition  is  required 
due  to  the  absence  of  upstream  influence.  Thus  the  boundary 
conditions  became: 


^ = ± Mx 

4>°=  4>x  E 0 


on  y 


as  x + -® 


(A.  9) 


In  each  of  the  above  cases,  the  gradient  boundary  condition 
is  treated  by  using  a one  sided  finite  difference  approxima- 


In  the  present  work,  the  farfield  boundary  conditions 
for  the  unsteady  perturbation  potential  are  treated  in  the 
same  manner  as  the  supersonic  steady  case.  The  wall  condition 
is  modified  to  allow  phase  changes  between  the  unsteady  per- 
turbation velocity  normal  to  the  wall  and  the  unsteady  pressure 
At  upstream  and  downstream  infinity  a zero  gradient  condition 
is  assumed  to  apply.  This  is  believed  to  be  accurate  for  the 
low  frequency  unsteady  conditions  considered,  and  the  adequacy 
of  the  approximation  has  in  fact  been  verified  for  the  quasi- 
steady (k=0)  case  for  which  a asymptotic  solution  exists 
(Equation  A. 7) . Thus  the  boundary  conditions  for  both  super- 
sonic and  subsonic  upstream  flows  are: 
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<j>  = + P e1^*  (<{>x  + i k4> ) on  y = ± h/2 

(A. 10) 

= 0 ,as  x ± » 


where  8*  is  the  phase  relationship  between  the  unsteady  velo- 
city normal  to  the  wall  and  the  unsteady  pressure  (-2  (<)>  +ik<)>) ) . 
It  is  recalled  that  the  ik<f>  term  in  the  above  expression  is 
reglected  in  the  low  frequency  approximation. 


A. 3 Three-Dimensional  Planar  Wing  in  Free  Air 


Figure  A.l  Three  Dimensional  Wing  Coordinate  System 


The  development  of  three  dimensional  subsonic  farfield 
approximations  for  the  steady  and  unsteady  perturbation  po- 
tentials proceeds  in  the  same  manner  as  described  above  for 
the  two  dimensional  flow.  As  before  the  method  involves  the 
approximation  of  various  integrals  over  the  winq  and  wake 
which  result  from  the  application  of  Green’s  theorem  to  the 
appropriate  partial  differential  equation.  Again,  Klunker 
has  used  the  method  to  develop  asymptotic  solutions  for  the 
three  dimensional  steady  flow  and  his  result  in  the  following 
form  is  used: 
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for  x + oo 


where  R = I (x2+K(y2+z2)  ] t(£,n)  is  the  wing  thickness  dis- 

tribution and  y(n)  is  the  spanwise  distribution  of  circulation. 

The  development  of  an  asymptotic  solution  for  the  un- 
steady perturbation  potential  follows  the  method  of  Klunker 
and  is  now  described  in  some  detail.  The  field  equation  for 
the  unsteady  perturbation  potential  (Equation  A. 10)  is  written 
as: 


L(<f>]  = K4>xx  +<l>yy  +<t>zz  “ 2ifi<J>x  + kflij) 


= {*x  *^x)  x 


(A. 12) 


The  application  of  Green's  theorem  to  the  linear  operator  L 
and  the  use  of  wing  and  wake  boundary  conditions  and  weak 
shock  conditions  results  in  the  following  integral  equation 
for  <t> : 


-151- 


(a.i: 


4>  (x,y , z) 


-JJ 

wiilg 


A<p(5»n)  ^ d£dn 


+f  Y (n>y  ip 5 e"lk(c"1)  d^dn 
span  1 

\ / 

wake  integral 


(4)|4>^)  ip;  d£dndr. 


where  ip  is  the  source  solution  to 
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The  use  of  the  source  function  ip  in  A.  13,  neglecting  the  v 
ume  integral  as  a higher  order  terra,  and  after  c°n®^e^l 
manipulation  and  approximation  (as  x +y  +z  •>»)  of  the 
integrals  results  in  the  following  farfield  solution. 
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Tt  in  the  equation,  is  an  integral  that  can  be  evaluated  using 
a rational  approximation  to  its  integrand  and  is  presented  at 
the  end  of  the  appendix. 

It  is  noted  that  both  the  steady  and  unsteady  farfield 
solutions  involve  integrals  over  the  wing  and  span  which  de- 
pend on  the  solution  (A<|>,y).  These  integrals  are  evaluated 
numerically  as  the  numerical  solution  proceeds  and  the  re- 
spective equations  are  used  to  update  the  values  of  the  steady 
or  unsteady  potential  on  the  farfield  boundaries. 

The  use  of  a rational  approximation  to  evaluate  the 
portion  of  the  wake  integral  given  as  Ij  above,  results  in 
the  following  function: 
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and  c = 0.372  with  defined  in  the  table 


V 

b 

V 

0 

1.0 

1 

-0.2418  6198 

2 

2.7918  027 

3 

-24.9910  79 

4 

111.5919  6 

5 

-271.4354  9 

6 

305.7528  8 

7 

| 41.1836  30 

8 

-545.9853  7 

9 

-644.7815  5 

10 

-328.7275  5 

11 

64.2795  11 

i 
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APPENDIX  B 

SUMMARY  TABLE  OF  UNSTEADY  AERODYNAMIC  FORCE  COEFFICIENTS 


mm. 


o +ae 


0 +BeJ 

J o 


Figure  B.l  Definition  of  Airfoil  Forces 

The  unsteady  aerodynamic  force  coefficients  calculated 
by  the  method  of  this  report  are  defined  and  tabulated  in  this 
appendix.  Results  are  presented  for  two  airfoils,  a NACA  64A00C 
and  NACA  64A410  undergoing  low  frequency  oscillations  (k*0+0.2) 
in  pitch,  plunge  and  control  surface  rotation  in  the  transonic 
speed  range  (M^®. 75+1. 2) . 

Figure  B.l  above  defines  the  airfoil  section  geometric 
and  force  parameters.  The  definition  of  aerodynamic  forces 
are  defined  positive  for  tail  down  pitch  and  control  surface 
oscillations  and  positive  for  vertical  translation  down.  The 
forces  are  given  by: 
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where  q is  freestream  dynamic  pressure  and  a,8,6l  are  the 
magnitu3es  of  the  unsteady  oscillations  in  pitch,  control 
surface  and  plunge.  Also,  w is  the  frequency  of  oscillation 
and  the  corresponding  reduced  frequency  is  defined  by  k=u>c/U 
It  is  noted  that  the  reduced  frequency  and  magnitude  of  plung- 
ing oscillation  are  based  on  airfoil  chord. 

Tables  of  unsteady  aerodynamic  coefficients,  so  defined, 
now  follow  for  k=0.0,  0.05,  0.1,  0.15,  0.2  for: 


• NACA  64A006  (pitch,  plunge,  control  surface  modes) 


= .8,  .85,  .9,  .95,  1.1,  1.15,  1.2 
“o  = °° 

Moo  = -85 

a„  = 0°,  1°,  2° 

• NACA  64A410  (pitch,  plunge) 


j = .75,  .8,  .85 

a0  = 0° 

M = .8 

a = 0°,  1°,  2° 
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HACA  64A006  (xu  = .75) 


ilACA  64A0O6  (xu  = .7D) (CONCLUDED) 
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